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

updating Z in QP roots

This commit is contained in:
Pierre-Francois Loos 2023-07-27 21:23:40 +02:00
parent 4fec6e88f5
commit 6b56acf66b
15 changed files with 50 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -15,5 +15,5 @@
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
T T F F F F
F F F T F F
# * unrestricted version available

2
input/options
View File

@ -9,7 +9,7 @@
# 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 F 0.0 F F
# ACFDT: AC Kx XBS

2
src/GF/G0F2.f90
View File

@ -78,7 +78,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF)
call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF,Z)
end if

13
src/GF/GF2_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF)
subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF,Z)
! Compute the graphical solution of the GF2 QP equation
@ -26,7 +26,7 @@ subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF)
! Output variables
double precision,intent(out) :: eGF(nBas)
double precision,intent(out) :: Z(nBas)
! Run Newton's algorithm to find the root
@ -48,22 +48,23 @@ subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF)
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
df = 1d0/(1d0 - dsigC)
w = w - f/df
w = w - df*f
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
end do
if(nIt == maxIt) then
write(*,*) 'Newton root search has not converged!'
eGF(p) = eHF(p)
else
eGF(p) = w
Z(p) = df
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGF(p)*HaToeV
write(*,*)

2
src/GF/evGF2.f90
View File

@ -101,7 +101,7 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF)
call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF,Z)
end if

2
src/GT/G0T0eh.f90
View File

@ -155,7 +155,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eHF,eGT)
call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eHF,eGT,Z)
end if

2
src/GT/G0T0pp.f90
View File

@ -164,7 +164,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
write(*,*)
call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,eHF,eGT)
Om1t,rho1t,Om2t,rho2t,eHF,eGT,Z)
end if

16
src/GT/GTeh_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
implicit none
include 'parameters.h'
@ -30,7 +30,9 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
double precision :: w
! Output variables
double precision,intent(out) :: eGT(nBas)
double precision,intent(out) :: Z(nBas)
sigC = 0d0
dsigC = 0d0
@ -54,20 +56,26 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
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
df = 1d0/(1d0 - dsigC)
w = w - f/df
w = w - df*f
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
end do
if(nIt == maxIt) then
write(*,*) 'Newton root search has not converged!'
eGT(p) = eGTlin(p)
else
eGT(p) = w
Z(p) = df
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV
write(*,*)
end if
end do

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(i,p,m,1)*rhoR(i,p,m,2)
num = rhoL(i,p,m,1)*rhoR(i,p,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
@ -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,2)*rhoR(p,a,m,1)
num = rhoL(p,a,m,1)*rhoR(p,a,m,1)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2

15
src/GT/GTpp_QP_graph.f90
View File

@ -1,5 +1,5 @@
subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,eGTlin,eGT)
Om1t,rho1t,Om2t,rho2t,eGTlin,eGT,Z)
implicit none
include 'parameters.h'
@ -33,7 +33,9 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
double precision :: w
! Output variables
double precision,intent(out) :: eGT(nBas)
double precision,intent(out) :: Z(nBas)
sigC = 0d0
dsigC = 0d0
@ -59,22 +61,29 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
dsigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
write (*,*) sigC
f = w - eHF(p) - sigC
df = 1d0 - dsigC
df = 1d0/(1d0 - dsigC)
w = w - f/df
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
end do
if(nIt == maxIt) then
eGT(p) = eGTlin(p)
write(*,*) 'Newton root search has not converged!'
else
eGT(p) = w
Z(p) = df
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV
write(*,*)
end if
end do
end subroutine

2
src/GT/evGTeh.f90
View File

@ -160,7 +160,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
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)
call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eOld,eGT,Z)
end if

2
src/GT/evGTpp.f90
View File

@ -191,7 +191,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
write(*,*)
call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,eOld,eGT)
Om1t,rho1t,Om2t,rho2t,eOld,eGT,Z)
end if

2
src/GW/G0W0.f90
View File

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

11
src/GW/GW_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW)
subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
! Compute the graphical solution of the QP equation
@ -34,6 +34,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW)
! Output variables
double precision,intent(out) :: eGW(nBas)
double precision,intent(out) :: Z(nBas)
! Run Newton's algorithm to find the root
@ -55,11 +56,11 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW)
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
df = 1d0/(1d0 - dsigC)
w = w - f/df
w = w - df*f
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
end do
@ -67,10 +68,12 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW)
if(nIt == maxIt) then
write(*,*) 'Newton root search has not converged!'
eGW(p) = eGWlin(p)
else
eGW(p) = w
Z(p) = df
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV
write(*,*)

2
src/GW/evGW.f90
View File

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