LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

fix memory leaks in ppBSE

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-27 10:11:35 +02:00
parent 91c6e7db40
commit 0a212628c3
16 changed files with 127 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -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

2
input/options
View File

@ -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

21
src/GF/G0F2.f90
View File

@ -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

27
src/GF/GF2_ppBSE2.f90
View File

@ -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)

12
src/GF/GF2_ppBSE2_dynamic_perturbation.f90
View File

@ -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)))

8
src/GF/GF2_self_energy.f90
View File

@ -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)

8
src/GF/GF2_self_energy_diag.f90
View File

@ -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)

19
src/GF/evGF2.f90
View File

@ -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

19
src/GF/qsGF2.f90
View File

@ -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

8
src/GF/regularized_self_energy_GF2.f90
View File

@ -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)

8
src/GF/regularized_self_energy_GF2_diag.f90
View File

@ -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)

33
src/GT/G0T0eh.f90
View File

@ -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)

34
src/GT/GTeh_excitation_density.f90
View File

@ -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

8
src/GT/GTeh_self_energy_diag.f90
View File

@ -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

3
src/GW/GW_ppBSE_static_kernel_C.f90
View File

@ -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

3
src/GW/GW_ppBSE_static_kernel_D.f90
View File

@ -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)