LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

root search for GTeh

This commit is contained in:
Pierre-Francois Loos 2023-07-27 16:43:15 +02:00
parent 866ed94153
commit 044058ac66
17 changed files with 255 additions and 73 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
T F F F F
F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F F F F
F F F T F F
# * unrestricted version available

6
input/options
View File

@ -9,10 +9,10 @@
# 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
256 0.00001 T 5 T 0.0 F F
256 0.00001 T 5 F 0.0 F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 T 0.0 F F
256 0.00001 T 5 F 0.0 F F
# ACFDT: AC Kx XBS
F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
F F T T T
F F F F T

4
src/GT/G0T0eh.f90
View File

@ -155,9 +155,11 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
else
write(*,*) ' *** Root search not yet implemented in G0T0eh *** '
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
end if
! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies

2
src/GT/G0T0pp.f90
View File

@ -192,7 +192,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,eGTlin,eGT)
call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,eGTlin,eGT)
end if

75
src/GT/GTeh_QP_graph.f90 Normal file
View File

@ -0,0 +1,75 @@
subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
implicit none
include 'parameters.h'
! Iput 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
double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS,2)
double precision,intent(in) :: rhoR(nBas,nBas,nS,2)
double precision,intent(in) :: eGTlin(nBas)
! Local variables
integer :: p
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: GTeh_SigC,GTeh_dSigC
double precision :: sigC,dsigC
double precision :: f,df
double precision :: w
! Output variables
double precision,intent(out) :: eGT(nBas)
sigC = 0d0
dsigC = 0d0
! Run Newton's algorithm to find the root
do p=nC+1,nBas-nR
write(*,*) '-----------------'
write(*,'(A10,I3)') 'Orbital ',p
write(*,*) '-----------------'
w = eGTlin(p)
nIt = 0
f = 1d0
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
do while (abs(f) > thresh .and. nIt < maxIt)
nIt = nIt + 1
sigC = GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR)
dsigC = GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR)
f = w - eHF(p) - sigC
df = 1d0 - dsigC
w = w - f/df
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
end do
if(nIt == maxIt) then
write(*,*) 'Newton root search has not converged!'
else
eGT(p) = w
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV
write(*,*)
end if
end do
end subroutine

53
src/GT/GTeh_SigC.f90 Normal file
View File

@ -0,0 +1,53 @@
double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
! 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) :: 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) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS,2)
double precision,intent(in) :: rhoR(nBas,nBas,nS,2)
! Local variables
integer :: i,a,m
double precision :: num,eps
! Initialize
GTeh_SigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do m=1,nS
eps = w - e(i) + Om(m)
num = rhoL(i,p,m,1)*rhoR(i,p,m,2)
GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do m=1,nS
eps = w - e(a) - Om(m)
num = rhoL(p,a,m,1)*rhoR(p,a,m,2)
GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
enddo
enddo
end function

53
src/GT/GTeh_dSigC.f90 Normal file
View File

@ -0,0 +1,53 @@
double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
! 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) :: 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) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS,2)
double precision,intent(in) :: rhoR(nBas,nBas,nS,2)
! Local variables
integer :: i,a,m
double precision :: num,eps
! Initialize
GTeh_dSigC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do m=1,nS
eps = w - e(i) + Om(m)
num = rhoL(i,p,m,1)*rhoR(i,p,m,2)
GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do m=1,nS
eps = w - e(a) - Om(m)
num = rhoL(p,a,m,1)*rhoR(p,a,m,2)
GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
end function

6
src/GT/GTeh_excitation_density.f90
View File

@ -51,12 +51,12 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
do m=1,nS
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,2) = rhoL(p,q,m,2) + ERI(p,j,b,q)*Y(jb,m) + ERI(p,b,j,q)*X(jb,m)
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)
+ (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X(jb,m) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*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)
+ (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y(jb,m) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X(jb,m)
enddo
enddo

4
src/GT/GTeh_self_energy_diag.f90
View File

@ -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(p,i,m,2)*rhoR(i,p,m,2)
num = rhoL(i,p,m,1)*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(a,p,m,1)*rhoR(p,a,m,1)
num = rhoL(p,a,m,1)*rhoR(p,a,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

14
src/GT/QP_graph_GT.f90 → src/GT/GTpp_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,eGTlin,eGT)
subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eGT)
implicit none
include 'parameters.h'
@ -14,9 +14,9 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
double precision,intent(in) :: eGTlin(nBas)
@ -26,7 +26,7 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: SigmaC_GT,dSigmaC_GT
double precision,external :: GTpp_SigC,GTpp_dSigC
double precision :: sigC,dsigC
double precision :: f,df
double precision :: w
@ -54,8 +54,8 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
nIt = nIt + 1
sigC = SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2)
dsigC = dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2)
sigC = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2)
dsigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2)
write (*,*) sigC
f = w - eHF(p) - sigC
df = 1d0 - dsigC
@ -76,4 +76,4 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
end do
end subroutine QP_graph_GT
end subroutine

16
src/GT/SigmaC_GT.f90 → src/GT/GTpp_SigC.f90
View File

@ -1,4 +1,4 @@
double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2)
double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
! Compute diagonal of the correlation part of the self-energy
@ -17,9 +17,9 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
integer,intent(in) :: nR
integer,intent(in) :: nOO,nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
@ -29,7 +29,7 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
! Initialize
SigmaC_GT = 0d0
GTpp_SigC = 0d0
!----------------------------------------------
! Occupied part of the T-matrix self-energy
@ -37,8 +37,8 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
do i=nC+1,nO
do cd=1,nVV
eps = w + e(i) - Omega1(cd)
SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
eps = w + e(i) - Om1(cd)
GTpp_SigC = GTpp_SigC + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
enddo
enddo
@ -48,8 +48,8 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
do a=nO+1,nBas-nR
do kl=1,nOO
eps = w + e(a) - Omega2(kl)
SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
eps = w + e(a) - Om2(kl)
GTpp_SigC = GTpp_SigC + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
enddo
enddo

16
src/GT/dSigmaC_GT.f90 → src/GT/GTpp_dSigC.f90
View File

@ -1,4 +1,4 @@
double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2)
double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
! Compute diagonal of the correlation part of the self-energy
@ -18,9 +18,9 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
@ -30,7 +30,7 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
! Initialize
dSigmaC_GT = 0d0
GTpp_dSigC = 0d0
!----------------------------------------------
! Occupied part of the T-matrix self-energy
@ -39,8 +39,8 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
do i=nC+1,nO
do cd=1,nVV
eps = w + e(i) - Omega1(cd)
dSigmaC_GT = dSigmaC_GT - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = w + e(i) - Om1(cd)
GTpp_dSigC= GTpp_dSigC - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
@ -52,8 +52,8 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
do a=nO+1,nBas-nR
do kl=1,nOO
eps = w + e(a) - Omega2(kl)
dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = w + e(a) - Om2(kl)
GTpp_dSigC = GTpp_dSigC - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo

2
src/GW/G0W0.f90
View File

@ -149,7 +149,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize)
call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize)
end if

10
src/GW/QP_graph.f90 → src/GW/GW_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
! Compute the graphical solution of the QP equation
@ -27,7 +27,7 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: SigmaC,dSigmaC
double precision,external :: GW_SigC,GW_dSigC
double precision :: sigC,dsigC
double precision :: f,df
double precision :: w
@ -53,8 +53,8 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
nIt = nIt + 1
sigC = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
dsigC = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
sigC = GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
dsigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
f = w - eHF(p) - SigC
df = 1d0 - dsigC
@ -80,4 +80,4 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
end do
end subroutine QP_graph
end subroutine

33
src/GW/SigmaC.f90 → src/GW/GW_SigC.f90
View File

@ -1,4 +1,4 @@
double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regularize)
double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
! Compute diagonal of the correlation part of the self-energy
@ -17,7 +17,7 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regular
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
logical,intent(in) :: regularize
@ -29,37 +29,36 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regular
! Initialize
SigmaC = 0d0
GW_SigC = 0d0
if (regularize) then
! Occupied part of the correlation self-energy
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
Dpijb = e(p) - e(i) + Omega(jb)
SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*(1d0-exp(-2d0*eta*Dpijb*Dpijb))/eps
eps = w - e(i) + Om(jb)
Dpijb = e(p) - e(i) + Om(jb)
GW_SigC = GW_SigC + 2d0*rho(p,i,jb)**2*(1d0-exp(-2d0*eta*Dpijb*Dpijb))/eps
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
Dpajb = e(p) - e(a) - Omega(jb)
SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(1d0-exp(-2d0*eta*Dpajb*Dpajb))/eps
eps = w - e(a) - Om(jb)
Dpajb = e(p) - e(a) - Om(jb)
GW_SigC = GW_SigC + 2d0*rho(p,a,jb)**2*(1d0-exp(-2d0*eta*Dpajb*Dpajb))/eps
enddo
enddo
! We add the static SRG term in the self-energy directly
! do i=nC+1,nO
! do jb=1,nS
! Dpijb = e(p) - e(i) + Omega(jb)
! Dpijb = e(p) - e(i) + Om(jb)
! SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*(exp(-2d0*eta*Dpijb*Dpijb)/Dpijb)
! enddo
! enddo
! do a=nO+1,nBas-nR
! do jb=1,nS
! Dpajb = e(p) - e(a) - Omega(jb)
! Dpajb = e(p) - e(a) - Om(jb)
! SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(exp(-2d0*eta*Dpajb*Dpajb)/Dpajb)
! enddo
! enddo
@ -68,17 +67,17 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regular
! Occupied part of the correlation self-energy
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
eps = w - e(i) + Om(jb)
GW_SigC = GW_SigC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
eps = w - e(a) - Om(jb)
GW_SigC = GW_SigC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
enddo
enddo
end if
end function SigmaC
end function

28
src/GW/dSigmaC.f90 → src/GW/GW_dSigC.f90
View File

@ -1,4 +1,4 @@
double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regularize)
double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
! Compute the derivative of the correlation part of the self-energy
@ -17,7 +17,7 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
logical,intent(in) :: regularize
@ -29,23 +29,23 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
! Initialize
dSigmaC = 0d0
GW_dSigC = 0d0
if (regularize) then
! Occupied part of the correlation self-energy
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
Dpijb = e(p) - e(i) + Omega(jb)
dSigmaC = dSigmaC - 2d0*rho(p,i,jb)**2*(1d0-exp(-2*eta*Dpijb*Dpijb))/(eps**2)
eps = w - e(i) + Om(jb)
Dpijb = e(p) - e(i) + Om(jb)
GW_dSigC = GW_dSigC - 2d0*rho(p,i,jb)**2*(1d0-exp(-2*eta*Dpijb*Dpijb))/(eps**2)
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
Dpajb = e(p) - e(a) - Omega(jb)
dSigmaC = dSigmaC - 2d0*rho(p,a,jb)**2*(1d0-exp(-2*eta*Dpajb*Dpajb))/(eps**2)
eps = w - e(a) - Om(jb)
Dpajb = e(p) - e(a) - Om(jb)
GW_dSigC = GW_dSigC - 2d0*rho(p,a,jb)**2*(1d0-exp(-2*eta*Dpajb*Dpajb))/(eps**2)
enddo
enddo
@ -56,8 +56,8 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = w - e(i) + Omega(jb)
dSigmaC = dSigmaC - 2d0*rho(p,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = w - e(i) + Om(jb)
GW_dSigC = GW_dSigC - 2d0*rho(p,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
enddo
@ -68,11 +68,11 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = w - e(a) - Omega(jb)
dSigmaC = dSigmaC - 2d0*rho(p,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = w - e(a) - Om(jb)
GW_dSigC = GW_dSigC - 2d0*rho(p,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
enddo
end if
end function dSigmaC
end function

2
src/GW/evGW.f90
View File

@ -161,7 +161,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize)
call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize)
end if