10
1
mirror of https://github.com/pfloos/quack synced 2025-05-06 07:05:33 +02:00

rm cG0W0 and cG0F2 method flags

This commit is contained in:
Loris Burth 2025-05-05 14:17:04 +02:00
parent c15d037fae
commit f6089556fe
20 changed files with 10 additions and 1301 deletions

View File

@ -18,8 +18,6 @@
F F F F
# G0T0eh evGTeh qsGTeh
F F F
# cG0W0 cG0F2
F F
# Parquet
F
# Rtest Utest Gtest

View File

@ -1,4 +1,4 @@
subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2, &
subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3, &
renorm,maxSCF, &
thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
@ -14,7 +14,6 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
logical,intent(in) :: dotest
logical,intent(in) :: doG0F2
logical,intent(in) :: docG0F2
logical,intent(in) :: doevGF2
logical,intent(in) :: doqsGF2
logical,intent(in) :: doufG0F02
@ -171,21 +170,4 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
end if
!------------------------------------------------------------------------
! Compute complex G0F2 electronic binding energies
!------------------------------------------------------------------------
if(docG0F2) then
call wall_time(start_GF)
call cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS, &
ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
end subroutine

View File

@ -1,152 +0,0 @@
subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF)
! Perform a one-shot second-order Green function calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: dophBSE
logical,intent(in) :: doppBSE
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nC
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: CAP(nOrb,nOrb)
! Local variables
integer :: p
double precision :: Ec
double precision :: EcBSE(nspin)
double precision,allocatable :: Re_SigC(:)
double precision,allocatable :: Im_SigC(:)
double precision,allocatable :: Re_Z(:)
double precision,allocatable :: Im_Z(:)
double precision,allocatable :: Re_eGFlin(:)
double precision, allocatable :: Im_eGFlin(:)
double precision,allocatable :: Re_eGF(:)
double precision,allocatable :: Im_eGF(:)
double precision, allocatable :: e_CAP(:)
! Hello world
write(*,*)
write(*,*)'*******************************'
write(*,*)'* Restricted G0F2 Calculation *'
write(*,*)'*******************************'
write(*,*)
! Memory allocation
allocate(Re_SigC(nOrb),Im_SigC(nOrb), Re_Z(nOrb),Im_Z(nOrb),&
Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),e_CAP(nOrb))
do p = 1, nOrb
e_CAP(p) = CAP(p,p)
end do
! Frequency-dependent second-order contribution
if(regularize) then
write(*,*) "Regularisation not implemented (yet)"
!call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
else
call cRGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_CAP)
end if
Re_eGFlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
Im_eGFlin(:) = e_CAP(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
if(linearize) then
write(*,*) '*** Quasiparticle energies obtained by linearization ***'
Re_eGF(:) = Re_eGFlin(:)
Im_eGF(:) = Im_eGFlin(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call cRGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,eHF,e_cap,Re_eGF,Im_eGF,Re_Z,Im_Z)
end if
! Print results
! call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_cRG0F2(nOrb,nO,eHF,e_CAP,Re_SigC,Im_SigC,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec)
! Perform BSE@GF2 calculation
!
! if(dophBSE) then
!
! call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,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 + ERHF + sum(EcBSE)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Perform ppBSE@GF2 calculation
!
! if(doppBSE) then
!
! call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
!
! EcBSE(2) = 3d0*EcBSE(2)
!
! 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) =',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Testing zone
!
! if(dotest) then
!
! call dump_test_value('R','G0F2 correlation energy',Ec)
! call dump_test_value('R','G0F2 HOMO energy',eGF(nO))
! call dump_test_value('R','G0F2 LUMO energy',eGF(nO+1))
!
! end if
!
! deallocate(SigC, Z, eGFlin, eGF)
!
end subroutine

View File

@ -1,59 +0,0 @@
double precision function cRGF2_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps
double precision :: eta_tilde
double precision :: num
! Initialize
cRGF2_Im_SigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = Re_w + eHF(a) - eHF(i) - eHF(j)
eta_tilde = eta - Im_w + e_cap(i) -( e_cap(a) - e_cap(j))
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
cRGF2_Im_SigC = cRGF2_Im_SigC + num*eta_tilde/(eps**2 + eta_tilde**2)
end do
end do
end do
! Virtual part of the correlation self-energy
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = Re_w + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i)
cRGF2_Im_SigC = cRGF2_Im_SigC - num*eta_tilde/(eps**2 + eta_tilde**2)
end do
end do
end do
end function

View File

@ -1,56 +0,0 @@
double precision function cRGF2_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps
double precision :: eta_tilde
double precision :: num
! Initialize
cRGF2_Im_dSigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = Re_w + eHF(a) - eHF(i) - eHF(j)
eta_tilde = eta - Im_w + e_cap(i) -( e_cap(a) - e_cap(j))
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
cRGF2_Im_dSigC = cRGF2_Im_dSigC - 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2
end do
end do
end do
! Virtual part of the correlation self-energy
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = Re_w + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i)
cRGF2_Im_dSigC = cRGF2_Im_dSigC + 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2
end do
end do
end do
end function

View File

@ -1,97 +0,0 @@
subroutine cRGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,Re_eOld,Im_eold,Re_eGF,Im_eGF,Re_Z,Im_Z)
! Compute the graphical solution of the GF2 QP equation
implicit none
include 'parameters.h'
! Input variables
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
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: Re_eGFlin(nBas)
double precision,intent(in) :: Im_eGFlin(nBas)
double precision,intent(in) :: Re_eOld(nBas)
double precision,intent(in) :: Im_eOld(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: p
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: cRGF2_Re_SigC,cRGF2_Im_SigC,cRGF2_Re_dSigC,cRGF2_Im_dSigC
double precision :: Re_SigC,Im_SigC,Re_dSigC,Im_dSigC
double precision :: Re_f,Im_f,Re_df,Im_df
double precision :: Re_w,Im_w
! Output variables
double precision,intent(out) :: Re_eGF(nBas),Im_eGF(nBas)
double precision,intent(out) :: Re_Z(nBas),Im_Z(nBas)
! Run Newton's algorithm to find the root
write(*,*)'-----------------------------------------------------'
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','Re(e_GFlin) (eV)','Re(e_GF) (eV)','Re(Z)'
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','Im(e_GFlin) (eV)','Im(e_GF) (eV)','Im(Z)'
write(*,*)'-----------------------------------------------------'
do p=nC+1,nBas-nR
Re_w = Re_eGFlin(p)
Im_w = Im_eGFlin(p)
nIt = 0
Re_f = 1d0
Im_f = 0d0
do while (abs(cmplx(Re_f,Im_f,kind=8)) > thresh .and. nIt < maxIt)
nIt = nIt + 1
Re_SigC = cRGF2_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI)
Im_SigC = cRGF2_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI)
Re_dSigC = cRGF2_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI)
Im_dSigC = cRGF2_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI)
Re_f = Re_w - eHF(p) - Re_SigC
Im_f = Im_w - e_cap(p) - Im_SigC
Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
Im_df = Im_dSigC/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
Re_w = Re_w - Re_df*Re_f + Im_df*Im_f
Im_w = Im_w - Re_f*Im_df - Re_df*Im_f
end do
if(nIt == maxIt) then
Re_eGF(p) = Re_eGFlin(p)
Im_eGF(p) = Im_eGFlin(p)
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGFlin(p)*HaToeV,Re_eGF(p)*HaToeV,Re_Z(p),'Cvg Failed!'
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Im_eGFlin(p)*HaToeV,Im_eGF(p)*HaToeV,Im_Z(p),'Cvg Failed!'
else
Re_eGF(p) = Re_w
Im_eGF(p) = Im_w
Re_Z(p) = Re_df
Im_Z(p) = Im_df
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Re_eGFlin(p)*HaToeV,Re_eGF(p)*HaToeV,Re_Z(p)
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGFlin(p)*HaToeV,Im_eGF(p)*HaToeV,Im_Z(p)
write(*,*)'-----------------------------------------------------'
end if
end do
end subroutine

View File

@ -1,58 +0,0 @@
double precision function cRGF2_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps
double precision :: eta_tilde
double precision :: num
! Initialize
cRGF2_Re_SigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = Re_w + eHF(a) - eHF(i) - eHF(j)
eta_tilde = eta - Im_w + e_cap(i) -(e_cap(a) - e_cap(j))
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
cRGF2_Re_SigC = cRGF2_Re_SigC + num*eps/(eps**2 + eta_tilde**2)
end do
end do
end do
! Virtual part of the correlation self-energy
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = Re_w + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i)
cRGF2_Re_SigC = cRGF2_Re_SigC + num*eps/(eps**2 + eta_tilde**2)
end do
end do
end do
end function

View File

@ -1,59 +0,0 @@
double precision function cRGF2_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps
double precision :: eta_tilde
double precision :: num
! Initialize
cRGF2_Re_dSigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = Re_w + eHF(a) - eHF(i) - eHF(j)
eta_tilde = eta - Im_w + e_cap(i) - (e_cap(a) - e_cap(j))
cRGF2_Re_dSigC = cRGF2_Re_dSigC -&
(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
end do
end do
end do
! Virtual part of the correlation self-energy
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = Re_w + eHF(i) - eHF(a) - eHF(b)
eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i)
cRGF2_Re_dSigC = cRGF2_Re_dSigC -&
(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
end do
end do
end do
end function

View File

@ -1,87 +0,0 @@
subroutine cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap)
! Compute diagonal part of the GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
! Input variables
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
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
integer :: p
double precision :: eps
double precision :: eta_tilde
double precision :: num
double precision,allocatable :: Re_DS(:)
double precision,allocatable :: Im_DS(:)
! Output variables
double precision,intent(out) :: Re_SigC(nBas)
double precision,intent(out) :: Im_SigC(nBas)
double precision,intent(out) :: Re_Z(nBas)
double precision,intent(out) :: Im_Z(nBas)
! Initialize
allocate(Re_DS(nBas),Im_DS(nBas))
Re_SigC(:) = 0d0
Im_SigC(:) = 0d0
Re_DS(:) = 0d0
Im_DS(:) = 0d0
! Compute GF2 self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = e(p) + e(a) - e(i) - e(j)
eta_tilde = eta - e_cap(p) + e_cap(i) - (e_cap(a) - e_cap(j))
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2)
Im_SigC(p) = Im_SigC(p) + num*eta_tilde/(eps**2 + eta_tilde**2)
Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
Im_DS(p) = Im_DS(p) - 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do
end do
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = e(p) + e(i) - e(a) - e(b)
eta_tilde = eta + e_cap(p) - e_cap(a) - e_cap(b) + e_cap(i)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2)
Im_SigC(p) = Im_SigC(p) - num*eta_tilde/(eps**2 + eta_tilde**2)
Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do
end do
end do
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
deallocate(Re_DS,Im_DS)
end subroutine

View File

@ -1,4 +1,4 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh,max_diis,doACFDT, &
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW)
@ -17,7 +17,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh
logical,intent(in) :: doqsGW
logical,intent(in) :: doufG0W0
logical,intent(in) :: doufGW
logical,intent(in) :: docG0W0
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
@ -93,22 +92,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh
end if
!------------------------------------------------------------------------
! Perform cG0W0 calculation
!------------------------------------------------------------------------
if(docG0W0) then
call wall_time(start_GW)
call cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,CAP_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0W0 = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform evGW calculation
!------------------------------------------------------------------------

View File

@ -1,242 +0,0 @@
subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF)
! Perform G0W0 calculation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: doppBSE
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: doSRG
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: CAP(nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: eHF(nOrb)
! Local variables
logical :: print_W = .false.
logical :: plot_self = .false.
logical :: dRPA_W
integer :: isp_W
integer :: p
double precision :: flow
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcGM
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Re_SigC(:)
double precision,allocatable :: Im_SigC(:)
double precision,allocatable :: Re_Z(:)
double precision,allocatable :: Im_Z(:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: Re_eGWlin(:)
double precision, allocatable :: Im_eGWlin(:)
double precision,allocatable :: Re_eGW(:)
double precision,allocatable :: Im_eGW(:)
double precision, allocatable :: e_cap(:)
! Hello world
write(*,*)
write(*,*)'***************************************'
write(*,*)'* Restricted complex G0W0 Calculation *'
write(*,*)'***************************************'
write(*,*)
! Spin manifold and TDA for dynamical screening
isp_W = 1
dRPA_W = .true.
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamical screening!'
write(*,*)
end if
! SRG regularization
flow = 500d0
if(doSRG) then
! Not implemented
write(*,*) '*** SRG regularized G0W0 scheme ***'
write(*,*) '!!! No SRG with cRG0W0 !!!'
write(*,*)
end if
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), &
Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),e_cap(nOrb))
do p = 1, nOrb
e_cap(p) = CAP(p,p)
end do
!-------------------!
! Compute screening !
!-------------------!
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
!--------------------------!
! Compute spectral weights !
!--------------------------!
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
!------------------------!
! Compute GW self-energy !
!------------------------!
call cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
! Linearized or graphical solution?
Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
Im_eGWlin(:) = e_cap(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
Re_eGW(:) = Re_eGWlin(:)
Im_eGW(:) = Im_eGWlin(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call cRGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin,eHF,e_cap,Re_eGW,Im_eGW,Re_Z,Im_Z)
end if
! Plot self-energy, renormalization factor, and spectral function
if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
! Cumulant expansion
! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z)
! Compute the RPA correlation energy
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,Re_eGW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
!--------------!
! Dump results !
!--------------!
call print_cRG0W0(nOrb,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM,CAP)
!---------------------------!
! Perform phBSE calculation !
!---------------------------!
!
! if(dophBSE) then
!
! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
! nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem
!
! if(doACFDT) then
!
! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
! end if
!
!!---------------------------!
!! Perform ppBSE calculation !
!!---------------------------!
!
! if(doppBSE) then
!
! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Testing zone
!
! if(dotest) then
!
! call dump_test_value('R','G0W0 correlation energy',EcRPA)
! call dump_test_value('R','G0W0 HOMO energy',Re_eGW(nO))
! call dump_test_value('R','G0W0 LUMO energy',Re_eGW(nO+1))
!
! end if
!
end subroutine

View File

@ -1,56 +0,0 @@
double precision function cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
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) :: Re_e(nBas)
double precision,intent(in) :: Im_e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,m
double precision :: num,eps,eta_tilde
! Initialize
cRGW_Im_SigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do m=1,nS
eps = Re_w - Re_e(i) + Om(m)
eta_tilde = eta - Im_w + Im_e(i)
num = 2d0*rho(p,i,m)**2
cRGW_Im_SigC = cRGW_Im_SigC + num*eta_tilde/(eps**2 + eta_tilde**2)
end do
end do
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do m=1,nS
eps = Re_w - Re_e(a) - Om(m)
eta_tilde = eta + Im_w - Im_e(a)
num = 2d0*rho(p,a,m)**2
cRGW_Im_SigC =cRGW_Im_SigC - num*eta_tilde/(eps**2 + eta_tilde**2)
end do
end do
end function

View File

@ -1,56 +0,0 @@
double precision function cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
! Compute the derivative of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
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) :: Re_e(nBas)
double precision,intent(in) :: Im_e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,m
double precision :: num,eps,eta_tilde
! Initialize
cRGW_Im_dSigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do m=1,nS
eps = Re_w - Re_e(i) + Om(m)
eta_tilde = eta - Im_w + Im_e(i)
num = 2d0*rho(p,i,m)**2
cRGW_Im_dSigC = cRGW_Im_dSigC - 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2
end do
end do
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do m=1,nS
eps = Re_w - Re_e(a) - Om(m)
eta_tilde = eta + Im_w - Im_e(a)
num = 2d0*rho(p,a,m)**2
cRGW_Im_dSigC = cRGW_Im_dSigC + 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2
end do
end do
end function

View File

@ -1,106 +0,0 @@
subroutine cRGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin, &
Re_eOld,Im_eOld,Re_eGW,Im_eGW,Re_Z,Im_Z)
! Compute the graphical solution of the QP equation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
logical,intent(in) :: doSRG
double precision,intent(in) :: eta
double precision,intent(in) :: flow
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: Re_eGWlin(nBas)
double precision,intent(in) :: Im_eGWlin(nBas)
double precision,external :: cRGW_Re_SigC,cRGW_Re_dSigC
double precision,external :: cRGW_Im_SigC,cRGW_Im_dSigC
double precision,intent(in) :: Re_eOld(nBas)
double precision,intent(in) :: Im_eOld(nBas)
! Local variables
integer :: p
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision :: Re_SigC,Re_dSigC
double precision :: Im_SigC,Im_dSigC
double precision :: Re_f,Im_f,Re_df,Im_df
double precision :: Re_w
double precision :: Im_w
! Output variables
double precision,intent(out) :: Re_eGW(nBas),Im_eGW(nBas)
double precision,intent(out) :: Re_Z(nBas),Im_Z(nBas)
! Run Newton's algorithm to find the root
write(*,*)'-----------------------------------------------------'
write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Re(e_GWlin) (eV)','Re(e_GW (eV))','Re(Z)'
write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Im(e_GWlin) (eV)','Im(e_GW (eV))','Im(Z)'
write(*,*)'-----------------------------------------------------'
do p=nC+1,nBas-nR
Re_w = Re_eGWlin(p)
Im_w = Im_eGWlin(p)
nIt = 0
Re_f = 1d0
Im_f = 1d0
do while (sqrt(Re_f**2+Im_f**2) > thresh .and. nIt < maxIt)
nIt = nIt + 1
Re_SigC = cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
Im_SigC = cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
Re_dSigC = cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
Im_dSigC = cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
Re_f = Re_w - eHF(p) - Re_SigC
Im_f = Im_w - e_cap(p) - Im_SigC
Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
Im_df = Im_dSigC/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
Re_w = Re_w - Re_df*Re_f + Im_df*Im_f
Im_w = Im_w - Re_f*Im_df - Re_df*Im_f
end do
if(nIt == maxIt) then
Re_eGW(p) = Re_eGWlin(p)
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p),'Cvg Failed!'
else
Re_eGW(p) = Re_w
Im_eGW(p) = Im_w
Re_Z(p) = Re_df
Im_Z(p) = Im_df
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p)
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGWlin(p)*HaToeV,Im_eGW(p)*HaToeV,Im_Z(p)
end if
write(*,*)'-----------------------------------------------------'
end do
write(*,*)
end subroutine

View File

@ -1,57 +0,0 @@
double precision function cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
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) :: Re_e(nBas)
double precision,intent(in) :: Im_e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,m
double precision :: num,eps
double precision :: eta_tilde
! Initialize
cRGW_Re_SigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do m=1,nS
eps = Re_w - Re_e(i) + Om(m)
eta_tilde = eta - Im_w + Im_e(i)
num = 2d0*rho(p,i,m)**2
cRGW_Re_SigC = cRGW_Re_SigC + num*eps/(eps**2 + eta_tilde**2)
end do
end do
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do m=1,nS
eps = Re_w - Re_e(a) - Om(m)
eta_tilde = eta + Im_w - Im_e(a)
num = 2d0*rho(p,a,m)**2
cRGW_Re_SigC = cRGW_Re_SigC + num*eps/(eps**2 + eta_tilde**2)
end do
end do
end function

View File

@ -1,56 +0,0 @@
double precision function cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
! Compute the derivative of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
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) :: Re_e(nBas)
double precision,intent(in) :: Im_e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,m
double precision :: num,eps,eta_tilde
! Initialize
cRGW_Re_dSigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do m=1,nS
eps = Re_w - Re_e(i) + Om(m)
eta_tilde = eta - Im_w + Im_e(i)
num = 2d0*rho(p,i,m)**2
cRGW_Re_dSigC = cRGW_Re_dSigC - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
end do
end do
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do m=1,nS
eps = Re_w - Re_e(a) - Om(m)
eta_tilde = eta + Im_w - Im_e(a)
num = 2d0*rho(p,a,m)**2
cRGW_Re_dSigC = cRGW_Re_dSigC - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
end do
end do
end function

View File

@ -1,102 +0,0 @@
subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z,e_cap)
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: e_cap(nBas)
! Local variables
integer :: i,a,p,m
double precision :: num,eps
double precision :: eta_tilde
double precision,allocatable :: Re_DS(:)
double precision,allocatable :: Im_DS(:)
! Output variables
double precision,intent(out) :: Re_Sig(nBas)
double precision,intent(out) :: Im_Sig(nBas)
double precision,intent(out) :: Re_Z(nBas)
double precision,intent(out) :: Im_Z(nBas)
double precision,intent(out) :: EcGM
! Initialize
allocate(Re_DS(nBas),Im_DS(nBas))
Re_Sig(:) = 0d0
Im_Sig(:) = 0d0
Re_DS(:) = 0d0
Im_DS(:) = 0d0
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do m=1,nS
eps = e(p) - e(i) + Om(m)
eta_tilde = eta - e_cap(p) + e_cap(i)
num = 2d0*rho(p,i,m)**2
Re_Sig(p) = Re_Sig(p) + num*eps/(eps**2 + eta_tilde**2)
Im_Sig(p) = Im_Sig(p) + num*eta_tilde/(eps**2 + eta_tilde**2)
Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
Im_DS(p) = Im_DS(p) - 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do m=1,nS
eps = e(p) - e(a) - Om(m)
eta_tilde = eta + e_cap(p) - e_cap(a)
num = 2d0*rho(p,a,m)**2
Re_Sig(p) = Re_Sig(p) + num*eps/(eps**2 + eta_tilde**2)
Im_Sig(p) = Im_Sig(p) - num*eta_tilde/(eps**2 + eta_tilde**2)
Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do
end do
! Galitskii-Migdal correlation energy
! MAYBE MODIFY THIS TERM
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do m=1,nS
eps = e(a) - e(i) + Om(m)
num = 4d0*rho(a,i,m)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do
end do
end do
! Compute renormalization factor from derivative
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
deallocate(Re_DS,Im_DS)
end subroutine

View File

@ -137,12 +137,14 @@ program QuAcK
doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
doG0T0eh,doevGTeh,doqsGTeh, &
docG0W0,docG0F2, &
doParquet, &
doRtest,doUtest,doGtest)
doCAP = docG0W0 .or. docG0F2 .or. docRHF ! Add different cases if they need CAP
docG0W0 = docG0W0 .or. (doG0W0 .and. docRHF)
docG0F2 = docG0F2 .or. (doG0F2 .and. docRHF)
! Determine complex function calls
doCAP = docRHF
docG0W0 = doG0W0 .and. docRHF
docG0F2 = doG0F2 .and. docRHF
!--------------------------!
! Read options for methods !
!--------------------------!

View File

@ -386,7 +386,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
if(doGF .and. .not. docRHF) then
call wall_time(start_GF)
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,renorm_GF,maxSCF_GF, &
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF, &
thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,complex_ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
@ -424,7 +424,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
if(doGW .and. .not. docRHF) then
call wall_time(start_GW)
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF_GW,thresh_GW,max_diis_GW, &
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF_GW,thresh_GW,max_diis_GW, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, &
V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW)

View File

@ -10,7 +10,6 @@ subroutine read_methods(working_dir, &
doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
doG0T0eh,doevGTeh,doqsGTeh, &
docG0W0,docG0F2, &
doParquet, &
doRtest,doUtest,doGtest)
@ -34,7 +33,6 @@ subroutine read_methods(working_dir, &
logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(out) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
logical,intent(out) :: doG0T0eh,doevGTeh,doqsGTeh
logical,intent(out) :: docG0W0,docG0F2
logical,intent(out) :: doParquet
logical,intent(out) :: doRtest,doUtest,doGtest
@ -207,17 +205,6 @@ subroutine read_methods(working_dir, &
if(ans2 == 'T') doevGTeh = .true.
if(ans3 == 'T') doqsGTeh = .true.
! Read Complex methods
docG0W0 = .false.
docG0F2 = .false.
read(1,*)
read(1,*) ans1,ans2
if(ans1 == 'T') docG0W0 = .true.
if(ans2 == 'T') docG0F2 = .true.
! Read coupled channels methods
doParquet = .false.