LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:53 +01:00
Code Releases Activity

clean up UGW

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-30 22:01:44 +02:00
parent 76230b7990
commit 2d61eaa29f
21 changed files with 139 additions and 127 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -15,5 +15,5 @@
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F F F F F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F F T F F F F F F F
# * unrestricted version available # * unrestricted version available

2
input/options
View File

@ -11,7 +11,7 @@
# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg
256 0.00001 T 5 T 0.0 F F 256 0.00001 T 5 T 0.0 F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 T 0.1 T F 256 0.00001 T 5 F 0.0 T F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA # BSE: phBSE phBSE2 ppBSE dBSE dTDA

4
src/GT/G0T0eh.f90
View File

@ -1,6 +1,6 @@
subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & subroutine G0T0eh(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, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI,dipole_int,PHF,cHF,eHF) ERI,dipole_int,eHF)
! Perform ehG0T0 calculation ! Perform ehG0T0 calculation
@ -37,8 +37,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
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)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
! Local variables ! Local variables

40
src/GT/G0T0pp.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & subroutine G0T0pp(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,PHF,cHF,eHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -32,9 +32,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables ! Local variables
@ -101,9 +99,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI_MO,Cpp) call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI_MO,Dpp) call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
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))
@ -123,9 +121,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI_MO,Cpp) call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI_MO,Dpp) call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
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))
@ -139,10 +137,10 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
!---------------------------------------------- !----------------------------------------------
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,nC,nO,nV,nR,nOOs,nVVs,ERI,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,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
@ -179,9 +177,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
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))
@ -192,9 +190,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
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))
@ -211,7 +209,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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,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,eHF,eGT,EcBSE) ERI,dipole_int,eHF,eGT,EcBSE)
if(exchange_kernel) then if(exchange_kernel) then
@ -247,7 +245,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
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,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,eHF,eGT,EcBSE) Om2t,X2t,Y2t,rho1t,rho2t,ERI,eHF,eGT,EcBSE)
if(exchange_kernel) then if(exchange_kernel) then
@ -273,7 +271,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
call GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, & call GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,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,eHF,eGT,EcBSE) ERI,dipole_int,eHF,eGT,EcBSE)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

16
src/GT/GT.f90
View File

@ -86,10 +86,10 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(unrestricted) then if(unrestricted) then
call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) dipole_int_aa,dipole_int_bb,cHF,epsHF)
else else
call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,PHF,cHF,epsHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if end if
call wall_time(end_GT) call wall_time(end_GT)
@ -107,12 +107,12 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then if(unrestricted) then
call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
PHF,cHF,epsHF) cHF,epsHF)
else else
call evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call evGTpp(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,EHF,ERI,dipole_int,PHF,cHF,epsHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if end if
call wall_time(end_GT) call wall_time(end_GT)
@ -157,7 +157,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
print*,'Unrestricted version of G0T0eh not yet implemented! Sorry.' print*,'Unrestricted version of G0T0eh not yet implemented! Sorry.'
else else
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,PHF,cHF,epsHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call soG0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if end if
call wall_time(end_GT) call wall_time(end_GT)
@ -178,7 +180,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
print*,'Unrestricted version of evGTeh not yet implemented! Sorry.' print*,'Unrestricted version of evGTeh not yet implemented! Sorry.'
else else
call evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & call evGTeh(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,EHF,ERI,dipole_int,PHF,cHF,epsHF) singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if end if
call wall_time(end_GT) call wall_time(end_GT)

3
src/GT/UG0T0pp.f90
View File

@ -1,7 +1,7 @@
subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, &
nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -34,7 +34,6 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
double precision,intent(in) :: EUHF double precision,intent(in) :: EUHF
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)

21
src/GT/evGTeh.f90
View File

@ -1,6 +1,5 @@
subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & subroutine evGTeh(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,PHF, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
cHF,eHF)
! Perform self-consistent eigenvalue-only ehGT calculation ! Perform self-consistent eigenvalue-only ehGT calculation
@ -36,10 +35,8 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
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) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables ! Local variables
@ -122,14 +119,14 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Compute screening ! Compute screening
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,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,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
! Compute spectral weights ! Compute spectral weights
call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
@ -224,7 +221,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! 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,eGW,eGW,EcBSE) ! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
! if(exchange_kernel) then ! if(exchange_kernel) then
@ -258,7 +255,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! end if ! end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) ! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
! write(*,*) ! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------' ! write(*,*)'-------------------------------------------------------------------------------'
@ -275,7 +272,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! if(ppBSE) then ! if(ppBSE) then
! call Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE) ! call Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcppBSE)
! write(*,*) ! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------' ! write(*,*)'-------------------------------------------------------------------------------'
@ -297,7 +294,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! !
! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) ! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW) ! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW)
! call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI) ! call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! !
! call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE) ! call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE)

26
src/GT/evGTpp.f90
View File

@ -1,5 +1,5 @@
subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,PHF,cHF,eHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@ -33,10 +33,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -129,9 +127,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
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) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
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))
@ -148,9 +146,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
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) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
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))
@ -164,10 +162,10 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
!---------------------------------------------- !----------------------------------------------
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,nC,nO,nV,nR,nOOs,nVVs,ERI,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,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nooS,nVVt,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nooS,nVVt,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
@ -233,7 +231,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
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,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,eGT,eGT,EcBSE) ERI,dipole_int,eGT,eGT,EcBSE)
if(exchange_kernel) then if(exchange_kernel) then
@ -269,7 +267,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,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,eGT,eGT,EcBSE)
if(exchange_kernel) then if(exchange_kernel) then

3
src/GT/evUGTpp.f90
View File

@ -1,7 +1,7 @@
subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa, & eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa, &
ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -35,7 +35,6 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
double precision,intent(in) :: EUHF double precision,intent(in) :: EUHF
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)

4
src/GW/G0W0.f90
View File

@ -1,5 +1,5 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & subroutine G0W0(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,dipole_int,PHF,cHF,eHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform G0W0 calculation ! Perform G0W0 calculation
@ -36,8 +36,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
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)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
! Local variables ! Local variables

12
src/GW/GW.f90
View File

@ -85,10 +85,10 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(unrestricted) then if(unrestricted) then
call UG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call UG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) dipole_int_aa,dipole_int_bb,cHF,epsHF)
else else
call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,PHF,cHF,epsHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if end if
call wall_time(end_GW) call wall_time(end_GW)
@ -106,12 +106,12 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
call wall_time(start_GW) call wall_time(start_GW)
if(unrestricted) then if(unrestricted) then
call evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, &
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S, &
PHF,cHF,epsHF) ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,epsHF)
else else
call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & call evGW(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,EHF,ERI,dipole_int,PHF,cHF,epsHF) singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if end if
call wall_time(end_GW) call wall_time(end_GW)

1
src/GW/GW_QP_graph.f90
View File

@ -62,7 +62,6 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
end do end do
if(nIt == maxIt) then if(nIt == maxIt) then

8
src/GW/GW_SigC.f90
View File

@ -23,7 +23,7 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Local variables ! Local variables
integer :: i,a,m integer :: i,a,m
double precision :: eps double precision :: num,eps
! Initialize ! Initialize
@ -34,7 +34,8 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do i=nC+1,nO do i=nC+1,nO
do m=1,nS do m=1,nS
eps = w - e(i) + Om(m) eps = w - e(i) + Om(m)
GW_SigC = GW_SigC + 2d0*rho(p,i,m)**2*eps/(eps**2 + eta**2) num = 2d0*rho(p,i,m)**2
GW_SigC = GW_SigC + num*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
@ -43,7 +44,8 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do m=1,nS do m=1,nS
eps = w - e(a) - Om(m) eps = w - e(a) - Om(m)
GW_SigC = GW_SigC + 2d0*rho(p,a,m)**2*eps/(eps**2 + eta**2) num = 2d0*rho(p,a,m)**2
GW_SigC = GW_SigC + num*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo

8
src/GW/GW_dSigC.f90
View File

@ -23,7 +23,7 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Local variables ! Local variables
integer :: i,a,m integer :: i,a,m
double precision :: eps double precision :: num,eps
! Initialize ! Initialize
@ -34,7 +34,8 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do i=nC+1,nO do i=nC+1,nO
do m=1,nS do m=1,nS
eps = w - e(i) + Om(m) eps = w - e(i) + Om(m)
GW_dSigC = GW_dSigC - 2d0*rho(p,i,m)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 num = 2d0*rho(p,i,m)**2
GW_dSigC = GW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo
enddo enddo
@ -43,7 +44,8 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do m=1,nS do m=1,nS
eps = w - e(a) - Om(m) eps = w - e(a) - Om(m)
GW_dSigC = GW_dSigC - 2d0*rho(p,a,m)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 num = 2d0*rho(p,a,m)**2
GW_dSigC = GW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo
enddo enddo

7
src/GW/UG0W0.f90
View File

@ -1,6 +1,6 @@
subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform unrestricted G0W0 calculation ! Perform unrestricted G0W0 calculation
@ -40,7 +40,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
! Local variables ! Local variables
@ -153,8 +152,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Find graphical solution of the QP equation ! Find graphical solution of the QP equation
do is=1,nspin do is=1,nspin
call UGW_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), & call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), &
OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) OmRPA,rho_RPA(:,:,:,is),eHF(:,is),eGW(:,is),Z(:,is))
end do end do
end if end if

21
src/GW/UGW_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
! Compute the graphical solution of the QP equation ! Compute the graphical solution of the QP equation
@ -15,7 +15,7 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS,nspin) double precision,intent(in) :: rho(nBas,nBas,nS,nspin)
double precision,intent(in) :: eGWlin(nBas) double precision,intent(in) :: eGWlin(nBas)
@ -24,9 +24,9 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
integer :: p integer :: p
integer :: nIt integer :: nIt
integer,parameter :: maxIt = 10 integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6 double precision,parameter :: thresh = 1d-6
double precision,external :: USigmaC,dUSigmaC double precision,external :: UGW_SigC,UGW_dSigC
double precision :: sigC,dsigC double precision :: sigC,dsigC
double precision :: f,df double precision :: f,df
double precision :: w double precision :: w
@ -34,6 +34,7 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
! Output variables ! Output variables
double precision,intent(out) :: eGW(nBas) double precision,intent(out) :: eGW(nBas)
double precision,intent(out) :: Z(nBas)
! Run Newton's algorithm to find the root ! Run Newton's algorithm to find the root
@ -52,15 +53,14 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
nIt = nIt + 1 nIt = nIt + 1
sigC = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) sigC = UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
dsigC = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) dsigC = UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
f = w - eHF(p) - sigC f = w - eHF(p) - sigC
df = 1d0 - dsigC df = 1d0/(1d0 - dsigC)
w = w - f/df w = w - df*f
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
end do end do
@ -72,6 +72,7 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
else else
eGW(p) = w eGW(p) = w
Z(p) = df
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV
write(*,*) write(*,*)

22
src/GW/USigmaC.f90 → src/GW/UGW_SigC.f90
View File

@ -1,4 +1,4 @@
double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Compute diagonal of the correlation part of the self-energy ! Compute diagonal of the correlation part of the self-energy
@ -17,32 +17,36 @@ double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
integer :: i,a,jb integer :: i,a,jb
double precision :: eps double precision :: num,eps
! Initialize ! Initialize
USigmaC = 0d0 UGW_SigC = 0d0
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
do i=nC+1,nO do i=nC+1,nO
do jb=1,nS do jb=1,nS
eps = w - e(i) + Omega(jb) eps = w - e(i) + Om(jb)
USigmaC = uSigmaC + rho(p,i,jb)**2*eps/(eps**2 + eta**2) num = rho(p,i,jb)**2
UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2)
end do end do
end do end do
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do jb=1,nS do jb=1,nS
eps = w - e(a) - Omega(jb) eps = w - e(a) - Om(jb)
USigmaC = USigmaC + rho(p,a,jb)**2*eps/(eps**2 + eta**2) num = rho(p,a,jb)**2
UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2)
end do end do
end do end do
end function USigmaC end function

20
src/GW/dUSigmaC.f90 → src/GW/UGW_dSigC.f90
View File

@ -1,4 +1,4 @@
double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Compute the derivative of the correlation part of the self-energy ! Compute the derivative of the correlation part of the self-energy
@ -17,24 +17,25 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
integer :: i,a,jb integer :: i,a,jb
double precision :: eps double precision :: num,eps
! Initialize ! Initialize
dUSigmaC = 0d0 UGW_dSigC = 0d0
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
do i=nC+1,nO do i=nC+1,nO
do jb=1,nS do jb=1,nS
eps = w - e(i) + Omega(jb) eps = w - e(i) + Om(jb)
dUSigmaC = dUSigmaC + rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2 num = rho(p,i,jb)**2
UGW_dSigC = UGW_dSigC + num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
@ -42,9 +43,10 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do jb=1,nS do jb=1,nS
eps = w - e(a) - Omega(jb) eps = w - e(a) - Om(jb)
dUSigmaC = dUSigmaC + rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 num = rho(p,a,jb)**2
UGW_dSigC = UGW_dSigC + num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end function dUSigmaC end function

17
src/GW/evGW.f90
View File

@ -1,6 +1,5 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & subroutine evGW(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,dipole_int,PHF, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
cHF,eHF)
! Perform self-consistent eigenvalue-only GW calculation ! Perform self-consistent eigenvalue-only GW calculation
@ -36,9 +35,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
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) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(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)
@ -147,17 +144,17 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
if(linearize) then if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*) write(*,*)
eGW(:) = eHF(:) + SigC(:) eGW(:) = eHF(:) + SigC(:)
else else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*) write(*,*)
call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eGW,Z) call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eGW,Z)
end if end if

27
src/GW/evUGW.f90
View File

@ -1,6 +1,6 @@
subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA, &
dBSE,dTDA,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform self-consistent eigenvalue-only GW calculation ! Perform self-consistent eigenvalue-only GW calculation
@ -24,6 +24,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta double precision,intent(in) :: eta
logical,intent(in) :: regularize logical,intent(in) :: regularize
@ -34,7 +35,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
integer,intent(in) :: nR(nspin) integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin) integer,intent(in) :: nS(nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas,nBas)
@ -156,7 +156,24 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Solve the quasi-particle equation ! ! Solve the quasi-particle equation !
!-----------------------------------! !-----------------------------------!
eGW(:,:) = eHF(:,:) + SigC(:,:) if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGW(:,:) = eHF(:,:) + SigC(:,:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
do is=1,nspin
call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), &
OmRPA,rho_RPA(:,:,:,is),eOld(:,is),eGW(:,is),Z(:,is))
end do
end if
! Convergence criteria ! Convergence criteria

2
src/GW/self_energy_correlation_SRG.f90
View File

@ -104,4 +104,4 @@ subroutine self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,
! end do ! end do
! end do ! end do
end subroutine self_energy_correlation_SRG end subroutine