10
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00

introduce nBas_MOs in RGW

This commit is contained in:
Abdallah Ammar 2024-08-29 00:00:41 +02:00
parent 0366561ce3
commit c06871a0ff
5 changed files with 217 additions and 144 deletions

View File

@ -1,7 +1,10 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! ---
subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF, thresh, max_diis, doACFDT, &
exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, &
linearize, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Restricted GW module
@ -43,7 +46,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -51,18 +54,18 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs)
double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
double precision,intent(in) :: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)
double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart)
double precision,intent(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart)
! Local variables
@ -76,11 +79,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0W0 = ',t_GW,' seconds'
write(*,*)
end if
@ -93,11 +96,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
singlet,triplet,linearize,eta,regularize,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGW = ',t_GW,' seconds'
write(*,*)
end if
@ -109,13 +112,14 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doqsGW) then
call wall_time(start_GW)
call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, &
ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, &
dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
write(*,*)
end if
@ -127,13 +131,15 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doSRGqsGW) then
call wall_time(start_GW)
call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, &
dophBSE, dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, &
nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, &
ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, &
PHF, cHF, eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
write(*,*)
end if
@ -145,11 +151,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufG0W0) then
call wall_time(start_GW)
call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
! TODO
call ufG0W0(dotest,TDA_W,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0W0 = ',t_GW,' seconds'
write(*,*)
end if
@ -161,11 +168,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufGW) then
call wall_time(start_GW)
call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
! TODO
call ufGW(dotest,TDA_W,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufGW = ',t_GW,' seconds'
write(*,*)
end if

View File

@ -1,6 +1,10 @@
subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA, &
dBSE,dTDA,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, &
BSE, BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nNuc, &
ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, &
X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -32,30 +36,30 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs)
double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
double precision,intent(inout):: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)
double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart)
double precision,intent(inout):: dipole_int_MO(nBas_MOs,nBas_MOs,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_AOs_Sq
integer :: ispin
integer :: ixyz
integer :: n_diis
@ -114,7 +118,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Stuff
nBasSq = nBas*nBas
nBas_AOs_Sq = nBas_AOs*nBas_AOs
! TDA for W
@ -132,9 +136,32 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Memory allocation
allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS), &
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGW(nBas_MOs))
allocate(eOld(nBas_MOs))
allocate(Z(nBas_MOs))
allocate(c(nBas_AOs,nBas_MOs))
allocate(cp(nBas_MOs,nBas_MOs))
allocate(Fp(nBas_MOs,nBas_MOs))
allocate(SigC(nBas_MOs,nBas_MOs))
allocate(P(nBas_AOs,nBas_AOs))
allocate(F(nBas_AOs,nBas_AOs))
allocate(J(nBas_AOs,nBas_AOs))
allocate(K(nBas_AOs,nBas_AOs))
allocate(error(nBas_AOs,nBas_AOs))
allocate(SigCp(nBas_AOs,nBas_AOs))
allocate(Aph(nS,nS))
allocate(Bph(nS,nS))
allocate(Om(nS))
allocate(XpY(nS,nS))
allocate(XmY(nS,nS))
allocate(rho(nBas_MOs,nBas_MOs,nS))
allocate(error_diis(nBas_AOs_Sq,max_diis))
allocate(F_diis(nBas_AOs_Sq,max_diis))
! Initialization
@ -162,11 +189,11 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Buid Hartree matrix
call wall_time(t1)
call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J)
call Hartree_matrix_AO_basis(nBas_AOs,P,ERI_AO,J)
! Compute exchange part of the self-energy
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
call exchange_matrix_AO_basis(nBas_AOs,P,ERI_AO,K)
call wall_time(t2)
tt=tt+t2-t1
@ -174,11 +201,11 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call wall_time(tao1)
do ixyz=1,ncart
call AOtoMO(nBas,nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
do ixyz = 1, ncart
call AOtoMO(nBas_AOs, nBas_MOs, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO)
call wall_time(tao2)
@ -188,8 +215,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call wall_time(tlr1)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR_A(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -203,13 +230,13 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call wall_time(tex1)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
call GW_excitation_density(nBas_MOs,nC,nO,nR,nS,ERI_MO,XpY,rho)
call wall_time(tex2)
tex=tex+tex2-tex1
call wall_time(tsrg1)
call SRG_self_energy(flow,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call SRG_self_energy(flow,nBas_MOs,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call wall_time(tsrg2)
@ -218,7 +245,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Make correlation self-energy Hermitian and transform it back to AO basis
call wall_time(tmo1)
call MOtoAO(nBas,nBas,S,c,SigC,SigCp)
call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp)
call wall_time(tmo2)
tmo = tmo + tmo2 - tmo1
! Solve the quasi-particle equation
@ -234,18 +261,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,error_diis,F_diis,error,F)
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
call diagonalize_matrix(nBas_MOs, cp, eGW)
c = matmul(X, cp)
call AOtoMO(nBas,nBas,c,SigCp,SigC)
call AOtoMO(nBas_AOs, nBas_MOs, c, SigCp, SigC)
! Compute new density matrix in the AO basis
@ -262,19 +289,19 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
ET = trace_matrix(nBas_AOs,matmul(P,T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
EV = trace_matrix(nBas_AOs,matmul(P,V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EJ = 0.5d0*trace_matrix(nBas_AOs,matmul(P,J))
! Exchange energy
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
Ex = 0.25d0*trace_matrix(nBas_AOs,matmul(P,K))
! Total energy
@ -282,8 +309,9 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, &
SigC, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGW, dipole)
end do
!------------------------------------------------------------------------
@ -300,6 +328,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis)
stop
end if
@ -313,17 +343,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Cumulant expansion
call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z)
call RGWC(dotest,eta,nBas_MOs,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z)
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis)
! Perform BSE calculation
if(BSE) then
call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
call GW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, &
nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE)
if(exchange_kernel) then
@ -357,7 +388,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE)
call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, BSE, singlet, triplet, &
eta, nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

View File

@ -1,4 +1,8 @@
subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
! ---
subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, &
Z, ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole)
! Print useful information about qsRGW calculation
@ -7,7 +11,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
@ -19,11 +23,11 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
double precision,intent(in) :: EcRPA
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: eGW(nBas_MOs)
double precision,intent(in) :: c(nBas_AOs,nBas_MOs)
double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs)
double precision,intent(in) :: Z(nBas_MOs)
double precision,intent(in) :: EqsGW
double precision,intent(in) :: dipole(ncart)
@ -59,7 +63,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
'|','#','|','e_HF (eV)','|','Sig_GW (eV)','|','Z','|','e_GW (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
do p=1,nBas_MOs
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
end do
@ -110,13 +114,13 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Restricted qsGW orbital coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas_AOs, nBas_MOs, c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Restricted qsGW orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eGW)
call vecout(nBas_MOs, eGW)
write(*,*)
end if

View File

@ -1,6 +1,10 @@
subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, &
ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, &
ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, &
ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -34,30 +38,30 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs)
double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_AOs)
double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
double precision,intent(inout):: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)
double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart)
double precision,intent(inout):: dipole_int_MO(nBas_MOs,nBas_MOs,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_AOs_Sq
integer :: ispin
integer :: ixyz
integer :: n_diis
@ -112,7 +116,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Stuff
nBasSq = nBas*nBas
nBas_AOs_Sq = nBas_AOs*nBas_AOs
! TDA for W
@ -130,10 +134,31 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Memory allocation
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), &
Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGW(nBas_MOs))
allocate(Z(nBas_MOs))
allocate(c(nBas_AOs,nBas_MOs))
allocate(cp(nBas_MOs,nBas_MOs))
allocate(Fp(nBas_MOs,nBas_MOs))
allocate(SigC(nBas_MOs,nBas_MOs))
allocate(P(nBas_AOs,nBas_AOs))
allocate(F(nBas_AOs,nBas_AOs))
allocate(J(nBas_AOs,nBas_AOs))
allocate(K(nBas_AOs,nBas_AOs))
allocate(err(nBas_AOs,nBas_AOs))
allocate(SigCp(nBas_AOs,nBas_AOs))
allocate(Aph(nS,nS))
allocate(Bph(nS,nS))
allocate(Om(nS))
allocate(XpY(nS,nS))
allocate(XmY(nS,nS))
allocate(rho(nBas_MOs,nBas_MOs,nS))
allocate(err_diis(nBas_AOs_Sq,max_diis))
allocate(F_diis(nBas_AOs_Sq,max_diis))
! Initialization
@ -160,38 +185,38 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Build Hartree-exchange matrix
call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J)
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
call Hartree_matrix_AO_basis(nBas_AOs, P, ERI_AO, J)
call exchange_matrix_AO_basis(nBas_AOs, P, ERI_AO, K)
! AO to MO transformation of two-electron integrals
do ixyz=1,ncart
call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
do ixyz = 1, ncart
call AOtoMO(nBas_AOs, nBas_MOs, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO)
! Compute linear response
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR_A(ispin, dRPA, nBas_MOs, nC, nO, nV, nR, nS, 1d0, eGW, ERI_MO, Aph)
if(.not.TDA_W) call phLR_B(ispin, dRPA, nBas_MOs, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phLR(TDA_W, nS, Aph, Bph, EcRPA, Om, XpY, XmY)
if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om)
! Compute correlation part of the self-energy
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
call GW_excitation_density(nBas_MOs, nC, nO, nR, nS, ERI_MO, XpY, rho)
if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
if(regularize) call GW_regularization(nBas_MOs, nC, nO, nV, nR, nS, eGW, Om, rho)
call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call GW_self_energy(eta, nBas_MOs, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z)
! Make correlation self-energy Hermitian and transform it back to AO basis
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO(nBas,nBas,S,c,SigC,SigCp)
call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp)
! Solve the quasi-particle equation
@ -205,19 +230,19 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
ET = trace_matrix(nBas_AOs, matmul(P, T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
EV = trace_matrix(nBas_AOs, matmul(P, V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
EK = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K))
! Total energy
@ -228,26 +253,27 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,F)
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
call AOtoMO(nBas,nBas,c,SigCp,SigC)
call diagonalize_matrix(nBas_MOs, cp, eGW)
c = matmul(X, cp)
call AOtoMO(nBas_AOs, nBas_MOs, c, SigCp, SigC)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole)
call print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, &
ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole)
end do
!------------------------------------------------------------------------
@ -264,19 +290,21 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis)
stop
end if
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis)
! Perform BSE calculation
if(dophBSE) then
call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
call GW_phBSE(dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, &
nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE)
if(exchange_kernel) then
@ -310,7 +338,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE)
call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, dophBSE, singlet, triplet, &
eta, nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -327,7 +356,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(doppBSE) then
call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE)
call GW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, &
nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, eGW, EcBSE)
EcBSE(2) = 3d0*EcBSE(2)

View File

@ -309,11 +309,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGW) then
call wall_time(start_GW)
! TODO
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF_GW, thresh_GW, max_diis_GW, &
doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, &
lin_GW, eta_GW, reg_GW, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, &
V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW