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 RGT

Browse Source
This commit is contained in:
Abdallah Ammar 2024-08-29 00:47:12 +02:00
parent c06871a0ff
commit c25e934e8b
6 changed files with 229 additions and 161 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
src/GT/RGT.f90
View File

@ -1,7 +1,11 @@
subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,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 RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, &
maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
doppBSE, TDA_T, 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)
! T-matrix module ! T-matrix module
@ -44,7 +48,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
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
@ -52,18 +56,18 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
integer,intent(in) :: nS integer,intent(in) :: 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(in) :: ERI_MO(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
@ -78,11 +82,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT) call wall_time(start_GT)
call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,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_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0pp = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0pp = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -95,11 +99,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT) call wall_time(start_GT)
call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,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_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTpp = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTpp = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -111,13 +115,13 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
if(doqsGTpp) then if(doqsGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, dBSE, &
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, & dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, &
PHF,cHF,eHF) 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_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTpp = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTpp = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -129,11 +133,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
if(doufG0T0pp) then if(doufG0T0pp) then
call wall_time(start_GT) call wall_time(start_GT)
call ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call ufG0T0pp(dotest,TDA_T,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0T0pp = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0T0pp = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -146,11 +150,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT) call wall_time(start_GT)
call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,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_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0eh = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0eh = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -163,11 +167,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT) call wall_time(start_GT)
call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,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_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTeh = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTeh = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -179,13 +183,14 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
if(doqsGTeh) then if(doqsGTeh) then
call wall_time(start_GT) call wall_time(start_GT)
call qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, &
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, & dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, &
PHF,cHF,eHF) 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_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTeh = ',t_GT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTeh = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if

24
src/GT/print_qsRGTeh.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
! ---
subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, &
Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole)
! Print one-electron energies and other stuff for qsGTeh ! Print one-electron energies and other stuff for qsGTeh
@ -7,7 +11,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,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
@ -19,11 +23,11 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: EcRPA(nspin)
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) :: eGT(nBas) double precision,intent(in) :: eGT(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) :: EqsGT double precision,intent(in) :: EqsGT
double precision,intent(in) :: dipole(ncart) double precision,intent(in) :: dipole(ncart)
@ -58,7 +62,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
'|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|' '|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|'
write(*,*)'-------------------------------------------------------------------------------' 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)') & 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),'|',eGT(p)*HaToeV,'|' '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|'
end do end do
@ -109,13 +113,13 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTeh MO coefficients' write(*,'(A32)') ' qsGTeh MO coefficients'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c) call matout(nBas_AOs, nBas_MOs, c)
write(*,*) write(*,*)
end if end if
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTeh MO energies' write(*,'(A32)') ' qsGTeh MO energies'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eGT) call vecout(nBas_MOs, eGT)
write(*,*) write(*,*)
end if end if

24
src/GT/print_qsRGTpp.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
! ---
subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, Z, &
ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole)
! Print one-electron energies and other stuff for qsGT ! Print one-electron energies and other stuff for qsGT
@ -7,7 +11,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,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
@ -19,11 +23,11 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: EcRPA(nspin)
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) :: eGT(nBas) double precision,intent(in) :: eGT(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) :: EqsGT double precision,intent(in) :: EqsGT
double precision,intent(in) :: dipole(ncart) double precision,intent(in) :: dipole(ncart)
@ -58,7 +62,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
'|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|'
write(*,*)'-------------------------------------------------------------------------------' 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)') & 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),'|',eGT(p)*HaToeV,'|' '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|'
end do end do
@ -109,13 +113,13 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTpp MO coefficients' write(*,'(A32)') ' qsGTpp MO coefficients'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c) call matout(nBas_AOs, nBas_MOs, c)
write(*,*) write(*,*)
end if end if
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTpp MO energies' write(*,'(A32)') ' qsGTpp MO energies'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eGT) call vecout(nBas_MOs, eGT)
write(*,*) write(*,*)
end if end if

110
src/GT/qsRGTeh.f90
View File

@ -1,6 +1,10 @@
subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,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 qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, &
dophBSE2, TDA_T, 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 GTeh calculation ! Perform a quasiparticle self-consistent GTeh calculation
@ -33,31 +37,31 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
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
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nS integer,intent(in) :: 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
logical :: dRPA = .false. logical :: dRPA = .false.
integer :: nSCF integer :: nSCF
integer :: nBasSq integer :: nBas_AOs_Sq
integer :: ispin integer :: ispin
integer :: n_diis integer :: n_diis
double precision :: ET double precision :: ET
@ -113,7 +117,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Stuff ! Stuff
nBasSq = nBas*nBas nBas_AOs_Sq = nBas_AOs*nBas_AOs
! TDA for T ! TDA for T
@ -131,9 +135,29 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Memory allocation ! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(Aph(nS,nS), Bph(nS,nS), Om(nS), XpY(nS,nS), XmY(nS,nS))
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) allocate(eGT(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(Sig(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(Sigp(nBas_AOs,nBas_AOs))
allocate(err(nBas_AOs,nBas_AOs))
allocate(err_diis(nBas_AOs_Sq,max_diis), F_diis(nBas_AOs_Sq,max_diis))
allocate(rhoL(nBas_MOs,nBas_MOs,nS), rhoR(nBas_MOs,nBas_MOs,nS))
! Initialization ! Initialization
@ -161,20 +185,20 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! 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 linear response ! Compute linear response
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) call phLR_A(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -182,17 +206,17 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) call GTeh_excitation_density(nBas_MOs,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR)
if(regularize) call GTeh_regularization(nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) if(regularize) call GTeh_regularization(nBas_MOs,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) call GTeh_self_energy(eta,nBas_MOs,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z)
! Make correlation self-energy Hermitian and transform it back to AO basis ! Make correlation self-energy Hermitian and transform it back to AO basis
Sig = 0.5d0*(Sig + transpose(Sig)) Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas,nBas,S,c,Sig,Sigp) call MOtoAO(nBas_AOs, nBas_MOs, S, c, Sig, Sigp)
! Solve the quasi-particle equation ! Solve the quasi-particle equation
@ -207,7 +231,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
if(max_diis > 1) then if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis) 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 end if
@ -215,9 +239,8 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
Fp = matmul(transpose(X),matmul(F,X)) Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGT) call diagonalize_matrix(nBas_MOs, cp, eGT)
c = matmul(X,cp) c = matmul(X,cp)
Sigp = matmul(transpose(c),matmul(Sigp,c))
! Compute new density matrix in the AO basis ! Compute new density matrix in the AO basis
@ -234,19 +257,19 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! 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))
! Total energy ! Total energy
@ -254,8 +277,9 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! 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_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) call print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, Sig, &
Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole)
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -272,19 +296,21 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) write(*,*)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis)
stop stop
end if end if
! Deallocate memory ! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,err,err_diis,F_diis) deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis)
! Perform BSE calculation ! Perform BSE calculation
! if(BSE) then ! if(BSE) then
! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & ! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
! eGT,eGT,EcBSE) ! eGT,eGT,EcBSE)
! if(exchange_kernel) then ! if(exchange_kernel) then
@ -319,7 +345,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! end if ! end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) ! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
! write(*,*) ! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------' ! write(*,*)'-------------------------------------------------------------------------------'

143
src/GT/qsRGTpp.f90
View File

@ -1,6 +1,9 @@
subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,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 qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, 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 GT calculation ! Perform a quasiparticle self-consistent GT calculation
@ -31,25 +34,26 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
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
integer,intent(in) :: 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 :: iblock integer :: iblock
integer :: n_diis integer :: n_diis
@ -119,7 +123,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Stuff ! Stuff
nBasSq = nBas*nBas nBas_AOs_Sq = nBas_AOs*nBas_AOs
! TDA for T ! TDA for T
@ -137,16 +141,30 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Memory allocation ! Memory allocation
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(eGT(nBas_MOs))
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas), & allocate(eOld(nBas_MOs))
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) allocate(Z(nBas_MOs))
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & allocate(c(nBas_AOs,nBas_MOs))
Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & allocate(Fp(nBas_MOs,nBas_MOs))
Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & allocate(cp(nBas_MOs,nBas_MOs))
Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & allocate(Sig(nBas_MOs,nBas_MOs))
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
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(Sigp(nBas_AOs,nBas_AOs))
allocate(error_diis(nBas_AOs_Sq,max_diis))
allocate(F_diis(nBas_AOs_Sq,max_diis))
allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nBas_MOs,nBas_MOs,nVVs))
allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nBas_MOs,nBas_MOs,nOOs))
allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nBas_MOs,nBas_MOs,nVVt))
allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nBas_MOs,nBas_MOs,nOOt))
! Initialization ! Initialization
@ -174,15 +192,15 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! 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 linear response ! Compute linear response
@ -191,9 +209,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) call ppLR_C(iblock,nBas_MOs,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) call ppLR_D(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) if(.not.TDA_T) call ppLR_B(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
@ -204,9 +222,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) call ppLR_C(iblock,nBas_MOs,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) call ppLR_D(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) if(.not.TDA_T) call ppLR_B(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
@ -218,24 +236,24 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
iblock = 3 iblock = 3
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) call GTpp_excitation_density(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
iblock = 4 iblock = 4
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) call GTpp_excitation_density(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
if(regularize) then if(regularize) then
call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) call GTpp_regularization(eta,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) call GTpp_regularization(eta,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
end if end if
call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & call GTpp_self_energy(eta,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
! Make correlation self-energy Hermitian and transform it back to AO basis ! Make correlation self-energy Hermitian and transform it back to AO basis
Sig = 0.5d0*(Sig + transpose(Sig)) Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas,nBas,S,c,Sig,Sigp) call MOtoAO(nBas_AOs, nBas_MOs, S, c, Sig, Sigp)
! Solve the quasi-particle equation ! Solve the quasi-particle equation
@ -249,22 +267,21 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
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,eGT) call diagonalize_matrix(nBas_MOs, cp, eGT)
c = matmul(X,cp) c = matmul(X, cp)
Sigp = matmul(transpose(c),matmul(Sigp,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
@ -277,19 +294,19 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! 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))
! Total energy ! Total energy
@ -297,8 +314,10 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! 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_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) call print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, &
eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, &
EqsGT, dipole)
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -315,19 +334,24 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) write(*,*)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis)
deallocate(Om1s, X1s, Y1s, rho1s)
deallocate(Om2s, X2s, Y2s, rho2s)
deallocate(Om1t, X1t, Y1t, rho1t)
deallocate(Om2t, X2t, Y2t, rho2t)
stop stop
end if end if
! Deallocate memory ! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis) deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis)
! Perform BSE calculation ! Perform BSE calculation
if(dophBSE) then if(dophBSE) then
call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas_MOs,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) ERI_MO,dipole_int_MO,eGT,eGT,EcBSE)
@ -363,7 +387,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
end if end if
call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas_MOs,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE) Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE)
@ -391,4 +415,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
end if end if
deallocate(Om1s, X1s, Y1s, rho1s)
deallocate(Om2s, X2s, Y2s, rho2s)
deallocate(Om1t, X1t, Y1t, rho1t)
deallocate(Om2t, X2t, Y2t, rho2t)
end subroutine end subroutine

10
src/QuAcK/RQuAcK.f90
View File

@ -330,11 +330,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGT) then if(doGT) then
call wall_time(start_GT) call wall_time(start_GT)
! TODO call RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, &
call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & maxSCF_GT, thresh_GT, max_diis_GT, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, &
maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & TDA_T, TDA, dBSE, dTDA, singlet, triplet, lin_GT, eta_GT, reg_GT, nNuc, ZNuc, rNuc, ENuc, &
TDA_T,TDA,dBSE,dTDA,singlet,triplet,lin_GT,eta_GT,reg_GT,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, &
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) dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT