ppBSE2@GF2

This commit is contained in:
Pierre-Francois Loos 2023-07-25 10:46:56 +02:00
parent 412e49f643
commit 91c6e7db40
12 changed files with 374 additions and 65 deletions

View File

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

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

View File

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

View File

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

View File

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

View File

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

164
src/GF/GF2_ppBSE2.f90 Normal file
View File

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

View File

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

View File

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

View File

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

View File

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

View File

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