PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-04 21:23:55 +01:00
Code Releases Activity

S^2 for ground state OK

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-10-03 22:16:38 +02:00
parent c15ea00351
commit ac6b92d169
13 changed files with 63 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.Be
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
1 2 2 0 0 1 3 1 0 0
# Znuc x y z # Znuc x y z
Be 0.0 0.0 0.0 Be 0.0 0.0 0.0

4
examples/molecule.H2
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
2 1 1 0 0 2 1 1 0 0
# Znuc x y z # Znuc x y z
H 0. 0. -0.7 H 0. 0. 0.0
H 0. 0. 0.7 H 0. 0. 1.4

2
examples/molecule.He
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
1 1 1 0 0 1 2 0 0 0
# Znuc x y z # Znuc x y z
He 0.0 0.0 0.0 He 0.0 0.0 0.0

4
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF MOM # RHF UHF MOM
T F F F T F
# MP2* MP3 MP2-F12 # MP2* MP3 MP2-F12
F F F F F F
# CCD CCSD CCSD(T) # CCD CCSD CCSD(T)
@ -7,7 +7,7 @@
# drCCD rCCD lCCD pCCD # drCCD rCCD lCCD pCCD
F F F F F F F F
# CIS* CIS(D) CID CISD # CIS* CIS(D) CID CISD
T T F F F F F F
# RPA* RPAx* ppRPA # RPA* RPAx* ppRPA
F F F F F F
# G0F2 evGF2 G0F3 evGF3 # G0F2 evGF2 G0F3 evGF3

6
input/options
View File

@ -5,14 +5,14 @@
# CC: maxSCF thresh DIIS n_diis # CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5 64 0.0000001 T 5
# spin: singlet triplet spin_conserved spin_flip TDA # spin: singlet triplet spin_conserved spin_flip TDA
T T T T T T T T F T
# GF: maxSCF thresh DIIS n_diis lin eta renorm # GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3 256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.0 F F F F F 256 0.00001 T 5 T 0.0 F F F F F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
T F T F F T
# BSE: BSE dBSE dTDA evDyn # BSE: BSE dBSE dTDA evDyn
T F T F F F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T 1000000 100000 10 0.3 10000 1234 T

2
src/CI/CIS.f90
View File

@ -55,6 +55,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
call diagonalize_matrix(nS,A,Omega) call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega) call print_excitation('CIS ',ispin,nS,Omega)
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A))
if(dump_trans) then if(dump_trans) then
print*,'Singlet CIS transition vectors' print*,'Singlet CIS transition vectors'
@ -82,6 +83,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
call diagonalize_matrix(nS,A,Omega) call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega) call print_excitation('CIS ',ispin,nS,Omega)
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A))
if(dump_trans) then if(dump_trans) then
print*,'Triplet CIS transition vectors' print*,'Triplet CIS transition vectors'

4
src/CI/UCIS.f90
View File

@ -80,6 +80,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
call diagonalize_matrix(nS_sc,A_sc,Omega_sc) call diagonalize_matrix(nS_sc,A_sc,Omega_sc)
call print_excitation('UCIS ',5,nS_sc,Omega_sc) call print_excitation('UCIS ',5,nS_sc,Omega_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
Omega_sc,transpose(A_sc),transpose(A_sc))
if(dump_trans) then if(dump_trans) then
print*,'Spin-conserved CIS transition vectors' print*,'Spin-conserved CIS transition vectors'
@ -122,6 +124,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
call quick_sort(Omega_sf,order(:),nS_sf) call quick_sort(Omega_sf,order(:),nS_sf)
call print_excitation('UCIS ',6,nS_sf,Omega_sf) call print_excitation('UCIS ',6,nS_sf,Omega_sf)
call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
Omega_sf,transpose(A_sf),transpose(A_sf))
if(dump_trans) then if(dump_trans) then
print*,'Spin-flip CIS transition vectors' print*,'Spin-flip CIS transition vectors'

2
src/HF/UHF.f90
View File

@ -233,6 +233,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH
! Compute final UHF energy ! Compute final UHF energy
call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
end subroutine UHF end subroutine UHF

16
src/HF/print_UHF.f90
View File

@ -1,4 +1,4 @@
subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
! Print one- and two-electron energies and other stuff for UHF calculation ! Print one- and two-electron energies and other stuff for UHF calculation
@ -7,6 +7,7 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin) integer,intent(in) :: nO(nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin) double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
@ -16,10 +17,14 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: Ex(nspin)
double precision,intent(in) :: EUHF double precision,intent(in) :: EUHF
integer :: i,j
integer :: ispin integer :: ispin
double precision :: HOMO(nspin) double precision :: HOMO(nspin)
double precision :: LUMO(nspin) double precision :: LUMO(nspin)
double precision :: Gap(nspin) double precision :: Gap(nspin)
double precision :: S2_exact
double precision :: S2
integer :: spin_state
! HOMO and LUMO ! HOMO and LUMO
@ -39,6 +44,11 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
end if end if
end do end do
S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(S,c(:,1:nO(2),2)))**2)
spin_state = nO(1) - nO(2) + 1
! Dump results ! Dump results
write(*,*) write(*,*)
@ -79,6 +89,10 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,I6)') ' 2S+1 :',spin_state
write(*,'(A40,F13.6)') ' <S**2> (exact) :',S2_exact
write(*,'(A40,F13.6)') ' <S**2> :',S2
write(*,'(A60)') '-------------------------------------------------'
write(*,*) write(*,*)
! Print results ! Print results

4
src/RPA/ACFDT.f90
View File

@ -94,7 +94,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
if(doXBS) then if(doXBS) then
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
@ -143,7 +143,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
if(doXBS) then if(doXBS) then
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)

4
src/RPA/print_transition_vectors.f90
View File

@ -108,8 +108,12 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,
end do end do
write(*,*) write(*,*)
print*,2d0*sum(X(:)**2 + Y(:)**2)
end do end do
! Thomas-Reiche-Kuhn sum rule
write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
write(*,*) write(*,*)

29
src/RPA/print_unrestricted_transition_vectors.f90
View File

@ -31,7 +31,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
integer :: ixyz integer :: ixyz
integer :: ispin integer :: ispin
integer,parameter :: maxS = 10 integer,parameter :: maxS = 10
double precision :: norm double precision :: S2
double precision,parameter :: thres_vec = 0.1d0 double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
@ -103,9 +103,20 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
print*,'---------------------------------------------' S2 = (nO(1) - nO(2))/2d0
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')' S2 = 2d0*S2+1d0
print*,'---------------------------------------------' S2 = 0.0d0
do jb=1,nSa
S2 = S2 + 4d0*(X(jb)**2 + Y(jb)**2)
end do
do jb=1,nSb
S2 = S2 - 4d0*(X(nSa+jb)**2 + Y(nSa+jb)**2)
end do
print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2
print*,'-------------------------------------------------------------'
! Spin-up transitions ! Spin-up transitions
@ -113,7 +124,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
do j=nC(1)+1,nO(1) do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
jb = jb + 1 jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)/sqrt(2d0) if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)
end do end do
end do end do
@ -121,7 +132,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
do j=nC(1)+1,nO(1) do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
jb = jb + 1 jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)/sqrt(2d0) if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)
end do end do
end do end do
@ -131,7 +142,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
jb = jb + 1 jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(jb)/sqrt(2d0) if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(nSa+jb)
end do end do
end do end do
@ -139,13 +150,15 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
jb = jb + 1 jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(jb)/sqrt(2d0) if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(nSa+jb)
end do end do
end do end do
write(*,*) write(*,*)
end do end do
! Thomas-Reiche-Kuhn sum rule
write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
write(*,*) write(*,*)

10
src/RPA/unrestricted_linear_response_A_matrix.f90
View File

@ -145,8 +145,8 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
jb = jb + 1 jb = jb + 1
! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B)' ! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B)'
A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) !& A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
! - (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i) - (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i)
end do end do
end do end do
@ -164,9 +164,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
jb = jb + 1 jb = jb + 1
A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) !& A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
! - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb) ! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb)
end do end do
end do end do