diff --git a/input/methods b/input/methods index 2305a6a..8848a96 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - F F F F T F + F F F F F F # * unrestricted version available diff --git a/input/options b/input/options index 316358c..7089880 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg 256 0.00001 T 5 T 0.0 F F # 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 F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 4c8ef31..b2d4900 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -1,6 +1,6 @@ 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, & - ERI,dipole_int,PHF,cHF,eHF) + ERI,dipole_int,eHF) ! 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) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) ! Local variables diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index f9ec259..9c382f1 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -1,5 +1,5 @@ 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) @@ -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) :: ERHF double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! 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)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI_MO,Dpp) + 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,Cpp) + 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)) @@ -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)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI_MO,Dpp) + 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,Cpp) + 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)) @@ -139,10 +137,10 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp !---------------------------------------------- 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 - 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, & 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)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + 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,Cpp) + 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)) @@ -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)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + 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,Cpp) + 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)) @@ -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, & 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 @@ -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, & 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 @@ -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, & 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(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GT/GT.f90 b/src/GT/GT.f90 index 702ebc8..eb48af6 100644 --- a/src/GT/GT.f90 +++ b/src/GT/GT.f90 @@ -86,10 +86,10 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(unrestricted) then 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, & - dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) + dipole_int_aa,dipole_int_bb,cHF,epsHF) else 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 call wall_time(end_GT) @@ -107,12 +107,12 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted call wall_time(start_GT) 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, & - PHF,cHF,epsHF) + cHF,epsHF) else 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 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.' else 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 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.' else 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 call wall_time(end_GT) diff --git a/src/GT/UG0T0pp.f90 b/src/GT/UG0T0pp.f90 index f5408fa..09c63cc 100644 --- a/src/GT/UG0T0pp.f90 +++ b/src/GT/UG0T0pp.f90 @@ -1,7 +1,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & spin_conserved,spin_flip,linearize,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) + dipole_int_aa,dipole_int_bb,cHF,eHF) ! 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) :: eHF(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_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 index f3a6792..c7f37bf 100644 --- a/src/GT/evGTeh.f90 +++ b/src/GT/evGTeh.f90 @@ -1,6 +1,5 @@ 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, & - cHF,eHF) + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! 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) :: nR integer,intent(in) :: nS - double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -122,14 +119,14 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! Compute screening - call phLR_A(ispin,dRPA,nBas,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) + 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,Bph) call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! 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 @@ -224,7 +221,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! 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 @@ -258,7 +255,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! 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(*,*)'-------------------------------------------------------------------------------' @@ -275,7 +272,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! 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(*,*)'-------------------------------------------------------------------------------' @@ -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,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) diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index 7ba89a8..d167f14 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -1,5 +1,5 @@ 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) @@ -33,10 +33,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) 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)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + 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,Cpp) + 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)) @@ -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)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + 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,Cpp) + 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)) @@ -164,10 +162,10 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T !---------------------------------------------- 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 - 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, & 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, & 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 @@ -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, & 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 diff --git a/src/GT/evUGTpp.f90 b/src/GT/evUGTpp.f90 index 2bb9c52..23b7e17 100644 --- a/src/GT/evUGTpp.f90 +++ b/src/GT/evUGTpp.f90 @@ -1,7 +1,7 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& 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) @@ -35,7 +35,6 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & double precision,intent(in) :: EUHF double precision,intent(in) :: eHF(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_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 6118758..8fec1fc 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,5 +1,5 @@ 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 @@ -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) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) ! Local variables diff --git a/src/GW/GW.f90 b/src/GW/GW.f90 index ff5692c..d2860bc 100644 --- a/src/GW/GW.f90 +++ b/src/GW/GW.f90 @@ -85,10 +85,10 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF if(unrestricted) then 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, & - dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) + dipole_int_aa,dipole_int_bb,cHF,epsHF) else 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 call wall_time(end_GW) @@ -106,12 +106,12 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF call wall_time(start_GW) if(unrestricted) then - call evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & - 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) + call evUGW(maxSCF,thresh,max_diis,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,dipole_int_aa,dipole_int_bb,cHF,epsHF) else 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 call wall_time(end_GW) diff --git a/src/GW/GW_QP_graph.f90 b/src/GW/GW_QP_graph.f90 index 0d1738e..d88b6a9 100644 --- a/src/GW/GW_QP_graph.f90 +++ b/src/GW/GW_QP_graph.f90 @@ -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 - end do if(nIt == maxIt) then diff --git a/src/GW/GW_SigC.f90 b/src/GW/GW_SigC.f90 index 0067c06..78564f4 100644 --- a/src/GW/GW_SigC.f90 +++ b/src/GW/GW_SigC.f90 @@ -23,7 +23,7 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Local variables integer :: i,a,m - double precision :: eps + double precision :: num,eps ! 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 m=1,nS 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 @@ -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 m=1,nS 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 diff --git a/src/GW/GW_dSigC.f90 b/src/GW/GW_dSigC.f90 index e31d9fc..19b14c7 100644 --- a/src/GW/GW_dSigC.f90 +++ b/src/GW/GW_dSigC.f90 @@ -23,7 +23,7 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Local variables integer :: i,a,m - double precision :: eps + double precision :: num,eps ! 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 m=1,nS 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 @@ -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 m=1,nS 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 diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 8f79754..8ab12d2 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -1,6 +1,6 @@ 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, & - dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) + dipole_int_aa,dipole_int_bb,cHF,eHF) ! 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) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) - double precision,intent(in) :: PHF(nBas,nBas,nspin) ! 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 do is=1,nspin - call UGW_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), & - OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) + call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), & + OmRPA,rho_RPA(:,:,:,is),eHF(:,is),eGW(:,is),Z(:,is)) end do end if diff --git a/src/GW/UGW_QP_graph.f90 b/src/GW/UGW_QP_graph.f90 index 1a5899c..5c5d9ca 100644 --- a/src/GW/UGW_QP_graph.f90 +++ b/src/GW/UGW_QP_graph.f90 @@ -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 @@ -15,7 +15,7 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) integer,intent(in) :: nS double precision,intent(in) :: eta 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) :: 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 :: nIt - integer,parameter :: maxIt = 10 + integer,parameter :: maxIt = 64 double precision,parameter :: thresh = 1d-6 - double precision,external :: USigmaC,dUSigmaC + double precision,external :: UGW_SigC,UGW_dSigC double precision :: sigC,dsigC double precision :: f,df 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 double precision,intent(out) :: eGW(nBas) + double precision,intent(out) :: Z(nBas) ! 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 - sigC = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) - dsigC = dUSigmaC(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 = UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho) 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 @@ -72,6 +72,7 @@ subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) else eGW(p) = w + Z(p) = df write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV write(*,*) diff --git a/src/GW/USigmaC.f90 b/src/GW/UGW_SigC.f90 similarity index 63% rename from src/GW/USigmaC.f90 rename to src/GW/UGW_SigC.f90 index c5e0817..1dc5d2d 100644 --- a/src/GW/USigmaC.f90 +++ b/src/GW/UGW_SigC.f90 @@ -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 @@ -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) :: nS 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) ! Local variables integer :: i,a,jb - double precision :: eps + double precision :: num,eps ! Initialize - USigmaC = 0d0 + UGW_SigC = 0d0 ! Occupied part of the correlation self-energy do i=nC+1,nO do jb=1,nS - eps = w - e(i) + Omega(jb) - USigmaC = uSigmaC + rho(p,i,jb)**2*eps/(eps**2 + eta**2) + eps = w - e(i) + Om(jb) + num = rho(p,i,jb)**2 + UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2) end do end do +! Virtual part of the correlation self-energy + do a=nO+1,nBas-nR do jb=1,nS - eps = w - e(a) - Omega(jb) - USigmaC = USigmaC + rho(p,a,jb)**2*eps/(eps**2 + eta**2) + eps = w - e(a) - Om(jb) + num = rho(p,a,jb)**2 + UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2) end do end do -end function USigmaC +end function diff --git a/src/GW/dUSigmaC.f90 b/src/GW/UGW_dSigC.f90 similarity index 65% rename from src/GW/dUSigmaC.f90 rename to src/GW/UGW_dSigC.f90 index 418e131..7e9ac65 100644 --- a/src/GW/dUSigmaC.f90 +++ b/src/GW/UGW_dSigC.f90 @@ -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 @@ -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) :: nS 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) ! Local variables integer :: i,a,jb - double precision :: eps + double precision :: num,eps ! Initialize - dUSigmaC = 0d0 + UGW_dSigC = 0d0 ! Occupied part of the correlation self-energy do i=nC+1,nO do jb=1,nS - eps = w - e(i) + Omega(jb) - dUSigmaC = dUSigmaC + rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2 + eps = w - e(i) + Om(jb) + num = rho(p,i,jb)**2 + UGW_dSigC = UGW_dSigC + num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 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 jb=1,nS - eps = w - e(a) - Omega(jb) - dUSigmaC = dUSigmaC + rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 + eps = w - e(a) - Om(jb) + 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 function dUSigmaC +end function diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 6cf308b..5fc6a09 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,6 +1,5 @@ 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, & - cHF,eHF) + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! 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) :: nR integer,intent(in) :: nS - double precision,intent(in) :: PHF(nBas,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) :: dipole_int(nBas,nBas,ncart) @@ -147,17 +144,17 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop if(linearize) then - write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' - write(*,*) + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) - eGW(:) = eHF(:) + SigC(:) + eGW(:) = eHF(:) + SigC(:) else - write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' - write(*,*) + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + 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 diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index a64d939..5cd8450 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -1,6 +1,6 @@ -subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & - dBSE,dTDA,spin_conserved,spin_flip,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) +subroutine evUGW(maxSCF,thresh,max_diis,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,dipole_int_aa,dipole_int_bb,cHF,eHF) ! 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) :: spin_conserved logical,intent(in) :: spin_flip + logical,intent(in) :: linearize double precision,intent(in) :: eta 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) :: nS(nspin) - double precision,intent(in) :: PHF(nBas,nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) 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 ! !-----------------------------------! - 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 diff --git a/src/GW/self_energy_correlation_SRG.f90 b/src/GW/self_energy_correlation_SRG.f90 index 2910269..deb00d1 100644 --- a/src/GW/self_energy_correlation_SRG.f90 +++ b/src/GW/self_energy_correlation_SRG.f90 @@ -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 subroutine self_energy_correlation_SRG +end subroutine