fix plot mistake

This commit is contained in:
Pierre-Francois Loos 2023-08-24 11:18:12 +02:00
parent 7efca56054
commit 9f7a9f5c21
6 changed files with 24 additions and 26 deletions

View File

@ -19,7 +19,7 @@ subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF,Z)
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: GF2_SigC,GF2_dSigC
double precision :: sigC,dsigC
double precision :: SigC,dSigC
double precision :: f,df
double precision :: w
@ -45,10 +45,10 @@ subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF,Z)
nIt = nIt + 1
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/(1d0 - dsigC)
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/(1d0 - dSigC)
w = w - df*f

View File

@ -149,7 +149,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
end if
call GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR)
call GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rhoL,rhoR)
! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies

View File

@ -29,7 +29,7 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: GTeh_SigC,GTeh_dSigC
double precision :: sigC,dsigC
double precision :: SigC,dSigC
double precision :: f,df
double precision :: w
@ -55,11 +55,10 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
nIt = nIt + 1
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/(1d0 - dsigC)
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/(1d0 - dSigC)
w = w - df*f
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f

View File

@ -28,7 +28,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: GTpp_SigC,GTpp_dSigC
double precision :: sigC,dsigC
double precision :: SigC,dSigC
double precision :: f,df
double precision :: w
@ -37,8 +37,8 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
double precision,intent(out) :: eGT(nBas)
double precision,intent(out) :: Z(nBas)
sigC = 0d0
dsigC = 0d0
SigC = 0d0
dSigC = 0d0
! Run Newton's algorithm to find the root
do p=nC+1,nBas-nR
@ -57,11 +57,10 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
nIt = nIt + 1
sigC = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
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/(1d0 - dsigC)
SigC = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
dSigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
f = w - eHF(p) - SigC
df = 1d0/(1d0 - dSigC)
w = w - f/df

View File

@ -140,7 +140,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
end if
call GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
call GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
! Compute the RPA correlation energy

View File

@ -28,7 +28,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: GW_SigC,GW_dSigC
double precision :: sigC,dsigC
double precision :: SigC,dSigC
double precision :: f,df
double precision :: w
@ -54,10 +54,10 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
nIt = nIt + 1
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/(1d0 - dsigC)
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/(1d0 - dSigC)
w = w - df*f