From 0a212628c3772e63a496193557ff708200fa0828 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 27 Jul 2023 10:11:35 +0200 Subject: [PATCH] fix memory leaks in ppBSE --- input/methods | 4 +-- input/options | 2 +- src/GF/G0F2.f90 | 21 ++++++++++--- src/GF/GF2_ppBSE2.f90 | 27 ++++++++-------- src/GF/GF2_ppBSE2_dynamic_perturbation.f90 | 12 ++++---- src/GF/GF2_self_energy.f90 | 8 +++-- src/GF/GF2_self_energy_diag.f90 | 8 +++-- src/GF/evGF2.f90 | 19 ++++++++++-- src/GF/qsGF2.f90 | 19 ++++++++++-- src/GF/regularized_self_energy_GF2.f90 | 8 +++-- src/GF/regularized_self_energy_GF2_diag.f90 | 8 +++-- src/GT/G0T0eh.f90 | 33 +++----------------- src/GT/GTeh_excitation_density.f90 | 34 +++++++++++++-------- src/GT/GTeh_self_energy_diag.f90 | 8 ++--- src/GW/GW_ppBSE_static_kernel_C.f90 | 3 +- src/GW/GW_ppBSE_static_kernel_D.f90 | 3 +- 16 files changed, 127 insertions(+), 90 deletions(-) diff --git a/input/methods b/input/methods index c89190f..da5c368 100644 --- a/input/methods +++ b/input/methods @@ -11,9 +11,9 @@ # phRPA* phRPAx* crRPA ppRPA F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 - F F F F F + T F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW - T F F F F F + F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh F F F F F F # * unrestricted version available diff --git a/input/options b/input/options index 0586628..174b10c 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 - T T F T T + F T T 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 diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index 2d94e4e..e32a9ca 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -62,11 +62,11 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu if(regularize) then - call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z) + call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,eHF,eHF,ERI,SigC,Z) else - call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z) + call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,eHF,ERI,SigC,Z) end if @@ -81,7 +81,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGFlin,ERI,eGF) + call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,eHF,eGFlin,ERI,eGF) end if @@ -110,6 +110,19 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu ! 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) + if(doppBSE) then + + call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if end subroutine diff --git a/src/GF/GF2_ppBSE2.f90 b/src/GF/GF2_ppBSE2.f90 index 910b935..a2ea9eb 100644 --- a/src/GF/GF2_ppBSE2.f90 +++ b/src/GF/GF2_ppBSE2.f90 @@ -1,4 +1,4 @@ -subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) +subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) ! Compute the Bethe-Salpeter excitation energies at the pp level @@ -19,7 +19,6 @@ subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, 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) @@ -80,9 +79,9 @@ subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! 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 GF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) + call GF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) + call GF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,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) @@ -100,9 +99,9 @@ subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! 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) + if(dBSE) & + call GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,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) @@ -132,9 +131,9 @@ subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! 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 GF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) + call GF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) + call GF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) @@ -153,9 +152,9 @@ subroutine GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! 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) + if(dBSE) & + call GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,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) diff --git a/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 b/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 index e3847d2..3c50d2a 100644 --- a/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 +++ b/src/GF/GF2_ppBSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eGF,ERI,dipole_int, & +subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,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 @@ -16,7 +16,6 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nO integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nS integer,intent(in) :: nOO integer,intent(in) :: nVV @@ -71,8 +70,9 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nO 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) +! if(.not.dTDA) call GF2_ppBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,OmBSE(ab),KB_dyn) + call GF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,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))) @@ -95,8 +95,8 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nO 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) +! if(.not.dTDA) call GF2_ppBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,OmBSE(ab),KB_dyn) + call GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,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))) diff --git a/src/GF/GF2_self_energy.f90 b/src/GF/GF2_self_energy.f90 index 106b2c6..3b98e66 100644 --- a/src/GF/GF2_self_energy.f90 +++ b/src/GF/GF2_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine GF2_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) +subroutine GF2_self_energy(eta,nBas,nC,nO,nV,nR,eHF,eGF2,ERI,SigC,Z) ! Compute GF2 self-energy and its renormalization factor @@ -8,7 +8,11 @@ subroutine GF2_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) ! Input variables double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF2(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) diff --git a/src/GF/GF2_self_energy_diag.f90 b/src/GF/GF2_self_energy_diag.f90 index b9f5aba..863babf 100644 --- a/src/GF/GF2_self_energy_diag.f90 +++ b/src/GF/GF2_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) +subroutine GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,eGF2,ERI,SigC,Z) ! Compute diagonal part of the GF2 self-energy and its renormalization factor @@ -8,7 +8,11 @@ subroutine GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) ! Input variables double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF2(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index e8deaf5..75feffa 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -82,11 +82,11 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr if(regularize) then - call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI,SigC,Z) + call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,eHF,eGF,ERI,SigC,Z) else - call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI,SigC,Z) + call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,eGF,ERI,SigC,Z) end if @@ -158,6 +158,19 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr ! 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) + if(doppBSE) then + + call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGF2 correlation energy (singlet) =',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGF2 correlation energy (triplet) =',3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGF2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGF2 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if end subroutine diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 2c1a1c1..6cc759a 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -149,11 +149,11 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr if(regularize) then - call regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI_MO,SigC,Z) + call regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,eHF,eGF,ERI_MO,SigC,Z) else - call GF2_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF,ERI_MO,SigC,Z) + call GF2_self_energy(eta,nBas,nC,nO,nV,nR,eHF,eGF,ERI_MO,SigC,Z) end if @@ -278,6 +278,19 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr ! 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) + if(doppBSE) then + + call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (singlet) =',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (triplet) =',3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 total energy =',ENuc + EqsGF2 + EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if end subroutine diff --git a/src/GF/regularized_self_energy_GF2.f90 b/src/GF/regularized_self_energy_GF2.f90 index c8df16d..df7146d 100644 --- a/src/GF/regularized_self_energy_GF2.f90 +++ b/src/GF/regularized_self_energy_GF2.f90 @@ -1,4 +1,4 @@ -subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) +subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,eHF,eGF2,ERI,SigC,Z) ! Compute GF2 self-energy and its renormalization factor @@ -8,7 +8,11 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC ! Input variables double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF2(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) diff --git a/src/GF/regularized_self_energy_GF2_diag.f90 b/src/GF/regularized_self_energy_GF2_diag.f90 index ce424cb..0e145cd 100644 --- a/src/GF/regularized_self_energy_GF2_diag.f90 +++ b/src/GF/regularized_self_energy_GF2_diag.f90 @@ -1,4 +1,4 @@ -subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) +subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,eHF,eGF2,ERI,SigC,Z) ! Compute diagonal part of the GF2 self-energy and its renormalization factor @@ -8,7 +8,11 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI ! Input variables double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF2(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 65a6c8c..10354cf 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -58,8 +58,8 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision,allocatable :: rhoL(:,:,:) - double precision,allocatable :: rhoR(:,:,:) + double precision,allocatable :: rhoL(:,:,:,:) + double precision,allocatable :: rhoR(:,:,:,:) double precision,allocatable :: eGT(:) double precision,allocatable :: eGTlin(:) @@ -102,42 +102,17 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),Sig(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & - rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),eGT(nBas),eGTlin(nBas)) + rhoL(nBas,nBas,nS,2),rhoR(nBas,nBas,nS,2),eGT(nBas),eGTlin(nBas)) !--------------------------------- -! Compute (singlet) RPA screening +! Compute (triplet) RPA screening !--------------------------------- -! allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),KA_sta(nS,nS),KB_sta(nS,nS)) - -! isp_W = 1 -! EcRPA = 0d0 - -! call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) -! if(.not.TDA_T) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - -! call phLR(TDA_T,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) -! call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - -! call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) -! call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) - -! deallocate(OmRPA,XpY_RPA,XmY_RPA,rho_RPA) - -!-------------------! -! Compute screening ! -!-------------------! - -! Spin manifold (triplet for GTeh) - ispin = 2 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) -! Aph(:,:) = Aph(:,:) + KA_sta(:,:) -! if(.not.TDA_T) Bph(:,:) = Bph(:,:) + KB_sta(:,:) - call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) if(print_T) call print_excitation('RPA@HF ',ispin,nS,Om) diff --git a/src/GT/GTeh_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90 index 3f3f867..553048f 100644 --- a/src/GT/GTeh_excitation_density.f90 +++ b/src/GT/GTeh_excitation_density.f90 @@ -14,17 +14,28 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) ! Local variables integer :: m,jb,p,q,j,b - double precision :: X,Y + double precision,allocatable :: X(:,:),Y(:,:) ! Output variables - double precision,intent(out) :: rhoL(nBas,nBas,nS) - double precision,intent(out) :: rhoR(nBas,nBas,nS) + double precision,intent(out) :: rhoL(nBas,nBas,nS,2) + double precision,intent(out) :: rhoR(nBas,nBas,nS,2) ! Initialization - rhoL(:,:,:) = 0d0 - rhoR(:,:,:) = 0d0 + rhoL(:,:,:,:) = 0d0 + rhoR(:,:,:,:) = 0d0 + + allocate(X(nS,nS),Y(nS,nS)) + + do m=1,nS + do jb=1,nS + X(jb,m) = 0.5d0*(XpY(m,jb) + XmY(m,jb)) + Y(jb,m) = 0.5d0*(XpY(m,jb) - XmY(m,jb)) + end do + end do + +! call matout(nS,nS,matmul(transpose(X),X) - matmul(transpose(Y),Y)) !$OMP PARALLEL & !$OMP SHARED(nC,nBas,nR,nO,nS,rhoL,rhoR,ERI,XpY,XmY) & @@ -39,14 +50,13 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) jb = jb + 1 do m=1,nS - X = 0.5d0*(XpY(m,jb) + XmY(m,jb)) - Y = 0.5d0*(XpY(m,jb) - XmY(m,jb)) + rhoL(p,q,m,1) = rhoL(p,q,m,1) + ERI(p,j,b,q)*X(jb,m) + ERI(p,b,j,q)*Y(jb,m) + rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(b,p,q,j)*Y(jb,m) + ERI(j,p,q,b)*X(jb,m) - rhoL(p,q,m) = rhoL(p,q,m) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y -! rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(p,b,j,q)*X + ERI(p,j,b,q)*Y - - rhoR(p,q,m) = rhoR(p,q,m) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y -! rhoR(p,q,m,2) = rhoR(p,q,m,2) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y + rhoR(p,q,m,1) = rhoR(p,q,m,1) & + + (2d0*ERI(b,q,p,j) - ERI(b,q,j,p))*X(jb,m) + (2d0*ERI(j,q,p,b) - ERI(j,q,b,p))*Y(jb,m) + rhoR(p,q,m,2) = rhoR(p,q,m,2) & + + (2d0*ERI(b,p,q,j) - ERI(b,p,j,q))*Y(jb,m) + (2d0*ERI(j,p,q,b) - ERI(j,p,b,q))*X(jb,m) enddo enddo diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90 index 1fcb612..58bb41e 100644 --- a/src/GT/GTeh_self_energy_diag.f90 +++ b/src/GT/GTeh_self_energy_diag.f90 @@ -16,8 +16,8 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rhoL(nBas,nBas,nS) - double precision,intent(in) :: rhoR(nBas,nBas,nS) + double precision,intent(in) :: rhoL(nBas,nBas,nS,2) + double precision,intent(in) :: rhoR(nBas,nBas,nS,2) ! Local variables @@ -46,7 +46,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do m=1,nS eps = e(p) - e(i) + Om(m) - num = rhoL(i,p,m)*rhoR(i,p,m) + num = rhoL(p,i,m,2)*rhoR(i,p,m,2) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -61,7 +61,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do m=1,nS eps = e(p) - e(a) - Om(m) - num = rhoL(p,a,m)*rhoR(p,a,m) + num = rhoL(a,p,m,1)*rhoR(p,a,m,1) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 diff --git a/src/GW/GW_ppBSE_static_kernel_C.f90 b/src/GW/GW_ppBSE_static_kernel_C.f90 index c3f94c1..bd015e3 100644 --- a/src/GW/GW_ppBSE_static_kernel_C.f90 +++ b/src/GW/GW_ppBSE_static_kernel_C.f90 @@ -1,4 +1,4 @@ -subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KC) +subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC) ! Compute the VVVV block of the static screening W for the pp-BSE @@ -14,7 +14,6 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: eta double precision,intent(in) :: lambda diff --git a/src/GW/GW_ppBSE_static_kernel_D.f90 b/src/GW/GW_ppBSE_static_kernel_D.f90 index 4948159..f15db35 100644 --- a/src/GW/GW_ppBSE_static_kernel_D.f90 +++ b/src/GW/GW_ppBSE_static_kernel_D.f90 @@ -1,4 +1,4 @@ -subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KD) +subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD) ! Compute the OOOO block of the static screening W for the pp-BSE @@ -15,7 +15,6 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda integer,intent(in) :: nR integer,intent(in) :: nS integer,intent(in) :: nOO - integer,intent(in) :: nVV double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)