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

introduce nBas_MOs in RGF

Browse Source
This commit is contained in:
Abdallah Ammar 2024-08-28 21:11:05 +02:00
parent d7403a946b
commit 0366561ce3
7 changed files with 144 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/GF/RG0F2.f90
View File

@ -51,7 +51,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
! Memory allocation ! Memory allocation
allocate(SigC(nBas),Z(nBas),eGFlin(nBas),eGF(nBas)) allocate(SigC(nBas), Z(nBas), eGFlin(nBas), eGF(nBas))
! Frequency-dependent second-order contribution ! Frequency-dependent second-order contribution
@ -133,4 +133,6 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
end if end if
deallocate(SigC, Z, eGFlin, eGF)
end subroutine end subroutine

70
src/GF/RGF.f90
View File

@ -1,7 +1,10 @@
subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & ! ---
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF) subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm, maxSCF, &
thresh, max_diis, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, linearize, &
eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, EHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF)
! Green's function module ! Green's function module
@ -39,7 +42,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
integer,intent(in) :: nBas integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nV integer,intent(in) :: nV
@ -47,18 +50,18 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: EHF double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas) double precision,intent(in) :: epsHF(nBas_MOs)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs)
double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs)
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart)
! Local variables ! Local variables
@ -71,12 +74,13 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doG0F2) then if(doG0F2) then
call wall_time(start_GF) call wall_time(start_GF)
call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & call RG0F2(dotest, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) linearize, eta, regularize, nBas_MOs, nC, nO, nV, nR, nS, &
ENuc, EHF, ERI_MO, dipole_int_MO, epsHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -89,12 +93,12 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
call wall_time(start_GF) call wall_time(start_GF)
call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & singlet,triplet,linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI,dipole_int,epsHF) ERI_MO,dipole_int_MO,epsHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -106,12 +110,14 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doqsGF2) then if(doqsGF2) then
call wall_time(start_GF) call wall_time(start_GF)
call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & call qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, &
nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, &
rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, EHF, S, &
X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGF2 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -123,11 +129,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doufG0F02) then if(doufG0F02) then
call wall_time(start_GF) call wall_time(start_GF)
call ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF) call ufRG0F02(dotest, nBas_MOs, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0F02 = ',t_GF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0F02 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -139,11 +145,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doG0F3) then if(doG0F3) then
call wall_time(start_GF) call wall_time(start_GF)
call RG0F3(dotest,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) call RG0F3(dotest, renorm, nBas_MOs, nC, nO, nV, nR, ERI_MO, epsHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -155,11 +161,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doevGF3) then if(doevGF3) then
call wall_time(start_GF) call wall_time(start_GF)
call evRGF3(dotest,maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nBas_MOs, nC, nO, nV, nR, ERI_MO, epsHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if

4
src/GF/evRGF2.f90
View File

@ -62,7 +62,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
! Memory allocation ! Memory allocation
allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) allocate(SigC(nBas), Z(nBas), eGF(nBas), eOld(nBas), error_diis(nBas,max_diis), e_diis(nBas,max_diis))
! Initialization ! Initialization
@ -189,4 +189,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
end if end if
deallocate(SigC, Z, eGF, eOld, error_diis, e_diis)
end subroutine end subroutine

24
src/GF/print_qsRGF2.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF,dipole)
! ---
subroutine print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, c, &
SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF, dipole)
! Print one-electron energies and other stuff for qsGF2 ! Print one-electron energies and other stuff for qsGF2
@ -7,17 +11,17 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,
! Input variables ! Input variables
integer,intent(in) :: nBas integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nSCF integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: Conv double precision,intent(in) :: Conv
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: eGF(nBas) double precision,intent(in) :: eGF(nBas_MOs)
double precision,intent(in) :: c(nBas) double precision,intent(in) :: c(nBas_AOs,nBas_MOs)
double precision,intent(in) :: SigC(nBas,nBas) double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs)
double precision,intent(in) :: Z(nBas) double precision,intent(in) :: Z(nBas_MOs)
double precision,intent(in) :: ET double precision,intent(in) :: ET
double precision,intent(in) :: EV double precision,intent(in) :: EV
double precision,intent(in) :: EJ double precision,intent(in) :: EJ
@ -53,7 +57,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,
'|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
do q=1,nBas do q = 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)') & 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)') &
'|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|' '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|'
end do end do
@ -102,12 +106,12 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGF2 MO coefficients' write(*,'(A32)') ' qsGF2 MO coefficients'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c) call matout(nBas_AOs, nBas_MOs, c)
write(*,*) write(*,*)
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGF2 MO energies' write(*,'(A32)') ' qsGF2 MO energies'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eGF) call matout(nBas_MOs, 1, eGF)
write(*,*) write(*,*)
end if end if

117
src/GF/qsRGF2.f90
View File

@ -1,6 +1,10 @@
subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,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 qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, &
dBSE, dTDA, 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 GF2 calculation ! Perform a quasiparticle self-consistent GF2 calculation
@ -29,25 +33,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas_AOs,nBas_MOs,nC,nO,nV,nR,nS
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs)
double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs)
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart)
double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart)
! Local variables ! Local variables
integer :: nSCF integer :: nSCF
integer :: nBasSq integer :: nBas_AOs_Sq
integer :: ispin integer :: ispin
integer :: n_diis integer :: n_diis
double precision :: EqsGF2 double precision :: EqsGF2
@ -94,7 +98,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Stuff ! Stuff
nBasSq = nBas*nBas nBas_AOs_Sq = nBas_AOs*nBas_AOs
! TDA ! TDA
@ -105,9 +109,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Memory allocation ! Memory allocation
allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(eGF(nBas_MOs))
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), & allocate(eOld(nBas_MOs))
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(c(nBas_AOs,nBas_MOs))
allocate(cp(nBas_MOs,nBas_MOs))
allocate(Fp(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(Z(nBas_MOs))
allocate(SigC(nBas_MOs,nBas_MOs))
allocate(SigCp(nBas_AOs,nBas_AOs))
allocate(error_diis(nBas_AOs_Sq,max_diis))
allocate(F_diis(nBas_AOs_Sq,max_diis))
! Initialization ! Initialization
@ -117,7 +139,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
Conv = 1d0 Conv = 1d0
P(:,:) = PHF(:,:) P(:,:) = PHF(:,:)
eOld(:) = eHF(:) eOld(:) = eHF(:)
eGF(:) = eHF(:) eGF(:) = eHF(:)
c(:,:) = cHF(:,:) c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0 F_diis(:,:) = 0d0
error_diis(:,:) = 0d0 error_diis(:,:) = 0d0
@ -135,25 +157,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Buid Hartree matrix ! Buid Hartree matrix
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 ! 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)
! AO to MO transformation of two-electron integrals ! AO to MO transformation of two-electron integrals
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 self-energy and renormalization factor ! Compute self-energy and renormalization factor
if(regularize) then if(regularize) then
call GF2_reg_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z) call GF2_reg_self_energy(eta, nBas_MOs, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
else else
call GF2_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z) call GF2_self_energy(eta, nBas_MOs, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
end if end if
@ -161,7 +183,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
SigC = 0.5d0*(SigC + transpose(SigC)) 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 ! Solve the quasi-particle equation
@ -169,28 +191,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) error = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
! DIIS extrapolation ! DIIS extrapolation
n_diis = min(n_diis+1,max_diis) n_diis = min(n_diis+1, max_diis)
if(abs(rcond) > 1d-7) then if(abs(rcond) > 1d-7) then
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)
else else
n_diis = 0 n_diis = 0
end if end if
! Diagonalize Hamiltonian in AO basis ! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X)) Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGF) call diagonalize_matrix(nBas_MOs, cp, eGF)
c = matmul(X,cp) c = matmul(X, cp)
SigCp = matmul(transpose(c),matmul(SigCp,c))
! Compute new density matrix in the AO basis ! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
! Save quasiparticles energy for next cycle ! Save quasiparticles energy for next cycle
@ -203,23 +224,23 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Kinetic energy ! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T)) ET = trace_matrix(nBas_AOs, matmul(P, T))
! Potential energy ! Potential energy
EV = trace_matrix(nBas,matmul(P,V)) EV = trace_matrix(nBas_AOs, matmul(P, V))
! Hartree energy ! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J))
! Exchange energy ! Exchange energy
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) Ex = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K))
! Correlation energy ! Correlation energy
call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec) call RMP2(.false., regularize, nBas_MOs, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec)
! Total energy ! Total energy
@ -230,8 +251,9 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Print results ! Print results
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole)
call print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) call print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, &
c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole)
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -248,19 +270,21 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
stop stop
end if end if
! Deallocate memory ! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis) deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
! Perform BSE calculation ! Perform BSE calculation
if(dophBSE) then if(dophBSE) then
call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE) call GF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, nC, nO, &
nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
@ -278,7 +302,8 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
if(doppBSE) then if(doppBSE) then
call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE) call GF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, &
nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

16
src/QuAcK/RQuAcK.f90
View File

@ -269,9 +269,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRPA) then if(doRPA) then
call wall_time(start_RPA) call wall_time(start_RPA)
! TODO call RRPA(dotest, dophRPA, dophRPAx, docrRPA, doppRPA, TDA, doACFDT, exchange_kernel, singlet, triplet, &
call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & nBas_MOs, nC, nO, nV, nR, nS, ENuc, ERHF, ERI_MO, dipole_int_MO, eHF)
nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S)
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
@ -289,11 +288,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGF) then if(doGF) then
call wall_time(start_GF) call wall_time(start_GF)
! TODO call RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm_GF, maxSCF_GF, &
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & thresh_GF, max_diis_GF, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, lin_GF, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, & eta_GF, reg_GF, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, &
nNuc,ZNuc,rNuc,ENuc,nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
@ -314,7 +312,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
! TODO ! TODO
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, & 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, & 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, & 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) ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GW) call wall_time(end_GW)

12
src/RPA/RRPA.f90
View File

@ -1,5 +1,5 @@
subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,cHF,S) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Random-phase approximation module ! Random-phase approximation module
@ -29,8 +29,6 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -49,7 +47,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds'
write(*,*) write(*,*)
end if end if
@ -65,7 +63,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPAx = ',t_RPA,' seconds'
write(*,*) write(*,*)
end if end if
@ -81,7 +79,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds'
write(*,*) write(*,*)
end if end if
@ -97,7 +95,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds'
write(*,*) write(*,*)
end if end if