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

fix root search

This commit is contained in:
Pierre-Francois Loos 2023-07-27 17:43:18 +02:00
parent 044058ac66
commit ee0cb0fd6a
18 changed files with 101 additions and 158 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/options
View File

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

12
src/GF/G0F2.f90
View File

@ -35,7 +35,6 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
double precision :: Ec
double precision :: EcBSE(nspin)
double precision,allocatable :: eGF(:)
double precision,allocatable :: eGFlin(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
@ -49,7 +48,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
! Memory allocation
allocate(SigC(nBas),Z(nBas),eGF(nBas),eGFlin(nBas))
allocate(SigC(nBas),Z(nBas),eGF(nBas))
if(linearize) then
@ -69,20 +68,17 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,eHF,ERI,SigC,Z)
end if
eGFlin(:) = eHF(:) + Z(:)*SigC(:)
if(linearize) then
eGF(:) = eGFlin(:)
eGF(:) = eHF(:) + Z(:)*SigC(:)
else
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,eHF,eGFlin,ERI,eGF)
call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF)
end if

17
src/GF/QP_graph_GF2.f90 → src/GF/GF2_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF)
! Compute the graphical solution of the GF2 QP equation
@ -10,7 +10,6 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF2lin(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
@ -19,14 +18,14 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: SigmaC_GF2,dSigmaC_GF2
double precision,external :: GF2_SigC,GF2_dSigC
double precision :: sigC,dsigC
double precision :: f,df
double precision :: w
! Output variables
double precision,intent(out) :: eGF2(nBas)
double precision,intent(out) :: eGF(nBas)
! Run Newton's algorithm to find the root
@ -37,7 +36,7 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
write(*,'(A10,I3)') 'Orbital ',p
write(*,*) '-----------------'
w = eGF2lin(p)
w = eHF(p)
nIt = 0
f = 1d0
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
@ -46,8 +45,8 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
nIt = nIt + 1
sigC = SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
dsigC = dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
sigC = GF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
dsigC = GF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
f = w - eHF(p) - sigC
df = 1d0 - dsigC
@ -64,9 +63,9 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
else
eGF2(p) = w
eGF(p) = w
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGF2(p)*HaToeV
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGF(p)*HaToeV
write(*,*)
end if

8
src/GF/SigmaC_GF2.f90 → src/GF/GF2_SigC.f90
View File

@ -1,4 +1,4 @@
double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
double precision function GF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
! Compute diagonal of the correlation part of the self-energy
@ -19,7 +19,7 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
integer :: i,j,a,b
double precision :: eps
SigmaC_GF2 = 0d0
GF2_SigC = 0d0
! Occupied part of the correlation self-energy
@ -28,7 +28,7 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
do a=nO+1,nBas-nR
eps = w + eHF(a) - eHF(i) - eHF(j)
SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
GF2_SigC = GF2_SigC + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
end do
end do
@ -41,7 +41,7 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
do b=nO+1,nBas-nR
eps = w + eHF(i) - eHF(a) - eHF(b)
SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
GF2_SigC = GF2_SigC + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
end do
end do

8
src/GF/dSigmaC_GF2.f90 → src/GF/GF2_dSigC.f90
View File

@ -1,4 +1,4 @@
double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
double precision function GF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
! Compute diagonal of the correlation part of the self-energy
@ -21,7 +21,7 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
! Initialize
dSigmaC_GF2 = 0d0
GF2_dSigC = 0d0
! Occupied part of the correlation self-energy
@ -30,7 +30,7 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
do a=nO+1,nBas-nR
eps = w + eHF(a) - eHF(i) - eHF(j)
dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
GF2_dSigC = GF2_dSigC - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
@ -43,7 +43,7 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
do b=nO+1,nBas-nR
eps = w + eHF(i) - eHF(a) - eHF(b)
dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
GF2_dSigC = GF2_dSigC - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do

17
src/GF/evGF2.f90
View File

@ -90,14 +90,19 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
end if
! Solve the quasi-particle equation
if(linearize) then
eGF(:) = eHF(:) + Z(:)*SigC(:)
else
eGF(:) = eHF(:) + SigC(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF)
end if
Conv = maxval(abs(eGF - eOld))

9
src/GT/G0T0eh.f90
View File

@ -62,7 +62,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
double precision,allocatable :: rhoR(:,:,:,:)
double precision,allocatable :: eGT(:)
double precision,allocatable :: eGTlin(:)
double precision,allocatable :: KA_sta(:,:)
double precision,allocatable :: KB_sta(:,:)
@ -102,7 +101,7 @@ 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,2),rhoR(nBas,nBas,nS,2),eGT(nBas),eGTlin(nBas))
rhoL(nBas,nBas,nS,2),rhoR(nBas,nBas,nS,2),eGT(nBas))
!---------------------------------
! Compute (triplet) RPA screening
@ -142,8 +141,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
! Solve the quasi-particle equation !
!-----------------------------------!
eGTlin(:) = eHF(:) + Z(:)*Sig(:)
! Linearized or graphical solution?
if(linearize) then
@ -151,14 +148,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGT(:) = eGTlin(:)
eGT(:) = eHF(:) + Z(:)*Sig(:)
else
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)
call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eHF,eGT)
end if

4
src/GT/GTeh_QP_graph.f90
View File

@ -51,8 +51,8 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
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)
sigC = GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGTlin,Om,rhoL,rhoR)
dsigC = GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGTlin,Om,rhoL,rhoR)
f = w - eHF(p) - sigC
df = 1d0 - dsigC

4
src/GT/GTeh_SigC.f90
View File

@ -35,7 +35,7 @@ double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
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)
num = rhoL(i,p,m,1)*rhoR(i,p,m,1)
GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
enddo
enddo
@ -45,7 +45,7 @@ double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
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)
num = rhoL(p,a,m,1)*rhoR(p,a,m,1)
GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
enddo
enddo

4
src/GT/GTeh_dSigC.f90
View File

@ -35,7 +35,7 @@ double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
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)
num = rhoL(i,p,m,1)*rhoR(i,p,m,1)
GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
@ -45,7 +45,7 @@ double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
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)
num = rhoL(p,a,m,1)*rhoR(p,a,m,1)
GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo

4
src/GT/GTeh_excitation_density.f90
View File

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

2
src/GT/GTeh_self_energy_diag.f90
View File

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

16
src/GT/evGTeh.f90
View File

@ -148,6 +148,22 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
eGT(:) = eHF(:) + Sig(:)
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGT(:) = eGT(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eOld,eGT)
end if
! Convergence criteria
Conv = maxval(abs(eGT - eOld))

12
src/GW/G0W0.f90
View File

@ -57,7 +57,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: eGWlin(:)
! Output variables
@ -94,8 +93,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
eGW(nBas),eGWlin(nBas))
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),eGW(nBas))
!-------------------!
! Compute screening !
@ -133,8 +131,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
! Solve the quasi-particle equation !
!-----------------------------------!
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
! Linearized or graphical solution?
if(linearize) then
@ -142,14 +138,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGW(:) = eGWlin(:)
eGW(:) = eHF(:) + Z(:)*SigC(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize)
call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eHF,eGW)
end if
@ -168,7 +164,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
! Deallocate memory
deallocate(SigC,Z,Om,XpY,XmY,rho,eGWlin)
deallocate(SigC,Z,Om,XpY,XmY,rho)
! Perform BSE calculation

9
src/GW/GW_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW)
! Compute the graphical solution of the QP equation
@ -15,11 +15,10 @@ subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regulari
integer,intent(in) :: nS
double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: eGWlin(nBas)
logical,intent(in) :: regularize
! Local variables
@ -53,8 +52,8 @@ subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regulari
nIt = nIt + 1
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)
sigC = GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
dsigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
f = w - eHF(p) - SigC
df = 1d0 - dsigC

65
src/GW/GW_SigC.f90
View File

@ -1,4 +1,4 @@
double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Compute diagonal of the correlation part of the self-energy
@ -19,65 +19,32 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regulariz
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
logical,intent(in) :: regularize
! Local variables
integer :: i,a,jb
integer :: i,a,m
double precision :: eps
double precision :: Dpijb,Dpajb
! Initialize
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) + 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) - 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
! Occupied part of the correlation self-energy
! 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) + 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) - Om(jb)
! SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(exp(-2d0*eta*Dpajb*Dpajb)/Dpajb)
! enddo
! enddo
else
! Occupied part of the correlation self-energy
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Om(jb)
GW_SigC = GW_SigC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
enddo
do i=nC+1,nO
do m=1,nS
eps = w - e(i) + Om(m)
GW_SigC = GW_SigC + 2d0*rho(p,i,m)**2*eps/(eps**2 + eta**2)
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Om(jb)
GW_SigC = GW_SigC + 2d0*rho(p,a,jb)**2*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)
GW_SigC = GW_SigC + 2d0*rho(p,a,m)**2*eps/(eps**2 + eta**2)
enddo
end if
enddo
end function

56
src/GW/GW_dSigC.f90
View File

@ -1,4 +1,4 @@
double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Compute the derivative of the correlation part of the self-energy
@ -19,60 +19,32 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regulari
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
logical,intent(in) :: regularize
! Local variables
integer :: i,j,a,b,jb
integer :: i,a,m
double precision :: eps
double precision :: Dpijb,Dpajb
! Initialize
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) + 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) - 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
! Occupied part of the correlation self-energy
else
! Occupied part of the correlation self-energy
do i=nC+1,nO
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
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
do i=nC+1,nO
do m=1,nS
eps = w - e(i) + Om(m)
GW_dSigC = GW_dSigC - 2d0*rho(p,i,m)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
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
do a=nO+1,nBas-nR
do m=1,nS
eps = w - e(a) - Om(m)
GW_dSigC = GW_dSigC - 2d0*rho(p,a,m)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
end if
enddo
end function

8
src/GW/evGW.f90
View File

@ -145,23 +145,19 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Solve the quasi-particle equation
eGW(:) = eHF(:) + SigC(:)
! Linearized or graphical solution?
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGW(:) = eGW(:)
eGW(:) = eHF(:) + SigC(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize)
call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eGW)
end if