From 91c6e7db40df54d82a5c953816c5421999591b6c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 25 Jul 2023 10:46:56 +0200 Subject: [PATCH] ppBSE2@GF2 --- input/methods | 4 +- input/options | 4 +- src/GF/G0F2.f90 | 38 +++-- src/GF/GF.f90 | 6 +- src/GF/GF2_phBSE2.f90 | 7 +- src/GF/GF2_phBSE2_dynamic_perturbation.f90 | 3 +- src/GF/GF2_ppBSE2.f90 | 164 +++++++++++++++++++++ src/GF/GF2_ppBSE2_dynamic_perturbation.f90 | 114 ++++++++++++++ src/GF/evGF2.f90 | 46 ++++-- src/GF/qsGF2.f90 | 42 +++--- src/GW/GW_ppBSE.f90 | 8 +- src/GW/GW_ppBSE_dynamic_perturbation.f90 | 3 +- 12 files changed, 374 insertions(+), 65 deletions(-) create mode 100644 src/GF/GF2_ppBSE2.f90 create mode 100644 src/GF/GF2_ppBSE2_dynamic_perturbation.f90 diff --git a/input/methods b/input/methods index 4ef876d..c89190f 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW - F F F F F F + T F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - F F F T F F + F F F F F F # * unrestricted version available diff --git a/input/options b/input/options index 0b1a49b..0586628 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + T T F T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 0 F # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg @@ -15,4 +15,4 @@ # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA - F F F F T + F F T T T diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index db1d272..2d94e4e 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -1,4 +1,4 @@ -subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & +subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a one-shot second-order Green function calculation @@ -8,7 +8,8 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & ! Input variables - logical,intent(in) :: BSE + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA @@ -33,8 +34,8 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & double precision :: Ec double precision :: EcBSE(nspin) - double precision,allocatable :: eGF2(:) - double precision,allocatable :: eGF2lin(:) + double precision,allocatable :: eGF(:) + double precision,allocatable :: eGFlin(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) @@ -48,7 +49,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGF2(nBas),eGF2lin(nBas)) + allocate(SigC(nBas),Z(nBas),eGF(nBas),eGFlin(nBas)) if(linearize) then @@ -69,33 +70,46 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & end if - eGF2lin(:) = eHF(:) + Z(:)*SigC(:) + eGFlin(:) = eHF(:) + Z(:)*SigC(:) if(linearize) then - eGF2(:) = eGF2lin(:) + eGF(:) = eGFlin(:) else write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2) + call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGFlin,ERI,eGF) end if ! Print results - call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF2,Ec) - call print_G0F2(nBas,nO,eHF,SigC,eGF2,Z,ENuc,ERHF,Ec) + call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) + call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec) ! Perform BSE2 calculation - if(BSE) then + if(dophBSE) then + + call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) - call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:)) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:)) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) end if +! Perform ppBSE2 calculation + + if(doppBSE) call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + end subroutine diff --git a/src/GF/GF.f90 b/src/GF/GF.f90 index 5297f27..5282856 100644 --- a/src/GF/GF.f90 +++ b/src/GF/GF.f90 @@ -81,7 +81,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,epsHF) else - call G0F2(dophBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & + call G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) end if call cpu_time(end_GF) @@ -104,7 +104,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,epsHF) else - call evGF2(dophBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & + call evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & ERI,dipole_int,epsHF) end if @@ -128,7 +128,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else - call qsGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & + call qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) end if call cpu_time(end_GF) diff --git a/src/GF/GF2_phBSE2.f90 b/src/GF/GF2_phBSE2.f90 index 0ebe1d3..dcce3cb 100644 --- a/src/GF/GF2_phBSE2.f90 +++ b/src/GF/GF2_phBSE2.f90 @@ -1,4 +1,4 @@ -subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE) +subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) ! Compute the second-order Bethe-Salpeter excitation energies @@ -20,7 +20,6 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -75,7 +74,7 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! Compute dynamic correction for BSE via perturbation theory if(dBSE) & - call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) + call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end if @@ -108,7 +107,7 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! Compute dynamic correction for BSE via perturbation theory if(dBSE) & - call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) + call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end if diff --git a/src/GF/GF2_phBSE2_dynamic_perturbation.f90 b/src/GF/GF2_phBSE2_dynamic_perturbation.f90 index dd91664..5bc5f66 100644 --- a/src/GF/GF2_phBSE2_dynamic_perturbation.f90 +++ b/src/GF/GF2_phBSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) +subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -19,7 +19,6 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER 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) :: eGF(nBas) double precision,intent(in) :: KA_sta(nS,nS) double precision,intent(in) :: KB_sta(nS,nS) diff --git a/src/GF/GF2_ppBSE2.f90 b/src/GF/GF2_ppBSE2.f90 new file mode 100644 index 0000000..910b935 --- /dev/null +++ b/src/GF/GF2_ppBSE2.f90 @@ -0,0 +1,164 @@ +subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + +! Compute the Bethe-Salpeter excitation energies at the pp level + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + + logical :: dRPA = .false. + + integer :: nOO + integer :: nVV + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + + double precision,allocatable :: Om1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + + double precision,allocatable :: Om2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + + double precision,allocatable :: KB_sta(:,:) + double precision,allocatable :: KC_sta(:,:) + double precision,allocatable :: KD_sta(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE(nspin) + +!------------------- +! Singlet manifold +!------------------- + + if(singlet) then + + write(*,*) '****************' + write(*,*) '*** Singlets ***' + write(*,*) '****************' + write(*,*) + + ispin = 1 + EcBSE(ispin) = 0d0 + + nOO = nO*(nO+1)/2 + nVV = nV*(nV+1)/2 + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & + KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + ! Compute BSE excitation energies + + if(.not.TDA) call GF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,KB_sta) + call GF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,KC_sta) + call GF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,KD_sta) + + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp) + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + + call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + + !----------------------------------------------------! + ! Compute the dynamical screening at the ppBSE level ! + !----------------------------------------------------! + +! if(dBSE) & +! call GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eGF,ERI,dipole_int, & +! Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) + + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + + end if + +!------------------- +! Triplet manifold +!------------------- + + if(triplet) then + + write(*,*) '****************' + write(*,*) '*** Triplets ***' + write(*,*) '****************' + write(*,*) + + ispin = 2 + EcBSE(ispin) = 0d0 + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & + KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + ! Compute BSE excitation energies + + if(.not.TDA) call GF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,KB_sta) + call GF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,KC_sta) + call GF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,KD_sta) + + + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp) + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + + call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + + !----------------------------------------------------! + ! Compute the dynamical screening at the ppBSE level ! + !----------------------------------------------------! + +! if(dBSE) & +! call GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eGF,ERI,dipole_int, & +! Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) + + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + + end if + +end subroutine diff --git a/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 b/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 new file mode 100644 index 0000000..e3847d2 --- /dev/null +++ b/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 @@ -0,0 +1,114 @@ +subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eGF,ERI,dipole_int, & + Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) + +! Compute dynamical effects via perturbation theory for BSE + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dTDA + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: Om2(nOO) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + double precision,intent(in) :: KB_sta(nVV,nOO) + double precision,intent(in) :: KC_sta(nVV,nVV) + double precision,intent(in) :: KD_sta(nOO,nOO) + +! Local variables + + integer :: ab,ij + + integer :: maxOO = 10 + integer :: maxVV = 10 + + double precision,allocatable :: Om1Dyn(:) + double precision,allocatable :: Om2Dyn(:) + double precision,allocatable :: Z1Dyn(:) + double precision,allocatable :: Z2Dyn(:) + + double precision,allocatable :: KB_dyn(:,:) + double precision,allocatable :: KC_dyn(:,:) + double precision,allocatable :: KD_dyn(:,:) + double precision,allocatable :: ZC_dyn(:,:) + double precision,allocatable :: ZD_dyn(:,:) + +! Memory allocation + + allocate(Om1Dyn(maxVV),Om2Dyn(maxOO),Z1Dyn(maxVV),Z2Dyn(maxOO), & + KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO), & + ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO)) + + if(dTDA) then + write(*,*) + write(*,*) '*** dynamical TDA activated ***' + write(*,*) + end if + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' First-order dynamical correction to static ppBSE2 double electron attachment energies ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + do ab=1,min(nVV,maxVV) + +! if(.not.dTDA) call GF2_ppBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGF,OmBSE(ab),KB_dyn) + call GF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGF,Om1(ab),KC_dyn,ZC_dyn) + + Z1Dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) + Om1Dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) + + Z1Dyn(ab) = 1d0/(1d0 - Z1Dyn(ab)) + Om1Dyn(ab) = Z1Dyn(ab)*Om1Dyn(ab) + + write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & + ab,Om1(ab)*HaToeV,(Om1(ab)+Om1Dyn(ab))*HaToeV,Om1Dyn(ab)*HaToeV,Z1Dyn(ab) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' First-order dynamical correction to static ppBSE2 double electron detachment energies ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + do ij=1,min(nOO,maxOO) + +! if(.not.dTDA) call GF2_ppBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGF,OmBSE(ab),KB_dyn) + call GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGF,Om2(ij),KD_dyn,ZD_dyn) + + Z2Dyn(ij) = dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij))) + Om2Dyn(ij) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) + + Z2Dyn(ij) = 1d0/(1d0 - Z2Dyn(ij)) + Om2Dyn(ij) = Z2Dyn(ij)*Om2Dyn(ij) + + write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & + ij,Om2(ij)*HaToeV,(Om2(ij)+Om2Dyn(ij))*HaToeV,Om2Dyn(ij)*HaToeV,Z2Dyn(ij) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + +end subroutine diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index d72e2e0..e8deaf5 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -1,4 +1,4 @@ -subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & +subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -8,7 +8,8 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & ! Input variables - logical,intent(in) :: BSE + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA @@ -41,7 +42,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & double precision :: EcBSE(nspin) double precision :: Conv double precision :: rcond - double precision,allocatable :: eGF2(:) + double precision,allocatable :: eGF(:) double precision,allocatable :: eOld(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) @@ -58,7 +59,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Initialization @@ -67,7 +68,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & n_diis = 0 e_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - eGF2(:) = eHF(:) + eGF(:) = eHF(:) eOld(:) = eHF(:) rcond = 0d0 @@ -81,39 +82,39 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & if(regularize) then - call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI,SigC,Z) else - call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI,SigC,Z) end if if(linearize) then - eGF2(:) = eHF(:) + Z(:)*SigC(:) + eGF(:) = eHF(:) + Z(:)*SigC(:) else - eGF2(:) = eHF(:) + SigC(:) + eGF(:) = eHF(:) + SigC(:) end if - Conv = maxval(abs(eGF2 - eOld)) + Conv = maxval(abs(eGF - eOld)) ! Print results - call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF2,Ec) - call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF2,ENuc,ERHF,Ec) + call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) + call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF-eOld,eGF) if(abs(rcond) < 1d-15) n_diis = 0 - eOld(:) = eGF2(:) + eOld(:) = eGF(:) ! Increment @@ -140,10 +141,23 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & ! Perform BSE2 calculation - if(BSE) then + if(dophBSE) then + + call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) - call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy =',sum(EcBSE(:)) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 total energy =',ENuc + EHF + sum(EcBSE(:)) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) end if +! Perform ppBSE2 calculation + + if(doppBSE) call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + end subroutine diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 8d27d52..2c1a1c1 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -1,4 +1,4 @@ -subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & +subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, & eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) @@ -12,7 +12,8 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & integer,intent(in) :: maxSCF integer,intent(in) :: max_diis double precision,intent(in) :: thresh - logical,intent(in) :: BSE + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA @@ -63,7 +64,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & double precision,allocatable :: F_diis(:,:) double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) - double precision,allocatable :: eGF2(:) + double precision,allocatable :: eGF(:) double precision,allocatable :: eOld(:) double precision,allocatable :: P(:,:) double precision,allocatable :: F(:,:) @@ -104,7 +105,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & ! Memory allocation - allocate(eGF2(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) @@ -116,7 +117,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & Conv = 1d0 P(:,:) = PHF(:,:) eOld(:) = eHF(:) - eGF2(:) = eHF(:) + eGF(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 @@ -148,11 +149,11 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & if(regularize) then - call regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) + call regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI_MO,SigC,Z) else - call GF2_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) + call GF2_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI_MO,SigC,Z) end if @@ -184,7 +185,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGF2) + call diagonalize_matrix(nBas,cp,eGF) c = matmul(X,cp) SigCp = matmul(transpose(c),matmul(SigCp,c)) @@ -194,8 +195,8 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & ! Save quasiparticles energy for next cycle - Conv = maxval(abs(eGF2 - eOld)) - eOld(:) = eGF2(:) + Conv = maxval(abs(eGF - eOld)) + eOld(:) = eGF(:) !------------------------------------------------------------------------ ! Compute total energy @@ -219,7 +220,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & ! Correlation energy - call MP2(regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF2,Ec) + call MP2(regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec) ! Total energy @@ -231,7 +232,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & !------------------------------------------------------------------------ call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,T,V,J,K,F,SigCp,Z, & + call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,P,T,V,J,K,F,SigCp,Z, & ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) enddo @@ -259,19 +260,24 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & ! Perform BSE calculation - if(BSE) then + if(dophBSE) then - call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGF2,EcBSE) + call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy =',sum(EcBSE(:)) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 total energy =',ENuc + EqsGF2 + sum(EcBSE(:)) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 correlation energy =',sum(EcBSE(:)) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 total energy =',ENuc + EqsGF2 + sum(EcBSE(:)) write(*,*)'-------------------------------------------------------------------------------' write(*,*) end if + +! Perform ppBSE2 calculation + + if(doppBSE) call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE) + end subroutine diff --git a/src/GW/GW_ppBSE.f90 b/src/GW/GW_ppBSE.f90 index 89d8f32..169dc4f 100644 --- a/src/GW/GW_ppBSE.f90 +++ b/src/GW/GW_ppBSE.f90 @@ -109,8 +109,8 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ! Compute BSE excitation energies if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KD_sta) + call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) @@ -161,8 +161,8 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ! Compute BSE excitation energies if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KD_sta) + call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) diff --git a/src/GW/GW_ppBSE_dynamic_perturbation.f90 b/src/GW/GW_ppBSE_dynamic_perturbation.f90 index 3befb3e..39afc87 100644 --- a/src/GW/GW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/GW_ppBSE_dynamic_perturbation.f90 @@ -56,7 +56,7 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO, ! Memory allocation - allocate(Om1Dyn(nVV),Om2Dyn(nOO),Z1Dyn(maxVV),Z2Dyn(maxOO), & + allocate(Om1Dyn(maxVV),Om2Dyn(maxOO),Z1Dyn(maxVV),Z2Dyn(maxOO), & KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO), & ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO)) @@ -77,7 +77,6 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO, ! if(.not.dTDA) call GW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,OmBSE(ab),KB_dyn) call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn) - Z1Dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) Om1Dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))