From 6b56acf66b0488beba051ad6165a4f0d2a0db71d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 27 Jul 2023 21:23:40 +0200 Subject: [PATCH] updating Z in QP roots --- input/methods | 2 +- input/options | 2 +- src/GF/G0F2.f90 | 2 +- src/GF/GF2_QP_graph.f90 | 13 +++++++------ src/GF/evGF2.f90 | 2 +- src/GT/G0T0eh.f90 | 2 +- src/GT/G0T0pp.f90 | 2 +- src/GT/GTeh_QP_graph.f90 | 16 ++++++++++++---- src/GT/GTeh_self_energy_diag.f90 | 4 ++-- src/GT/GTpp_QP_graph.f90 | 15 ++++++++++++--- src/GT/evGTeh.f90 | 2 +- src/GT/evGTpp.f90 | 2 +- src/GW/G0W0.f90 | 2 +- src/GW/GW_QP_graph.f90 | 11 +++++++---- src/GW/evGW.f90 | 2 +- 15 files changed, 50 insertions(+), 29 deletions(-) diff --git a/input/methods b/input/methods index bb2e24b..4ef876d 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index dd8f88a..f2119d9 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index f223287..f25a5d8 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -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 diff --git a/src/GF/GF2_QP_graph.f90 b/src/GF/GF2_QP_graph.f90 index 9899a46..aa39c77 100644 --- a/src/GF/GF2_QP_graph.f90 +++ b/src/GF/GF2_QP_graph.f90 @@ -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(*,*) diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index bf85965..776d3b3 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -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 diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 11b22f1..0c9d2b5 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -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 diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 92082ee..b5da52e 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -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 diff --git a/src/GT/GTeh_QP_graph.f90 b/src/GT/GTeh_QP_graph.f90 index f525298..86338af 100644 --- a/src/GT/GTeh_QP_graph.f90 +++ b/src/GT/GTeh_QP_graph.f90 @@ -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 diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90 index 753f5b3..6f3b913 100644 --- a/src/GT/GTeh_self_energy_diag.f90 +++ b/src/GT/GTeh_self_energy_diag.f90 @@ -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 diff --git a/src/GT/GTpp_QP_graph.f90 b/src/GT/GTpp_QP_graph.f90 index 1e0949b..ffa0368 100644 --- a/src/GT/GTpp_QP_graph.f90 +++ b/src/GT/GTpp_QP_graph.f90 @@ -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 diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 index 9e85485..f3a6792 100644 --- a/src/GT/evGTeh.f90 +++ b/src/GT/evGTeh.f90 @@ -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 diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index 1f1fdb6..7ba89a8 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -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 diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index ba3e7b1..f78febe 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -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 diff --git a/src/GW/GW_QP_graph.f90 b/src/GW/GW_QP_graph.f90 index d862284..0d1738e 100644 --- a/src/GW/GW_QP_graph.f90 +++ b/src/GW/GW_QP_graph.f90 @@ -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(*,*) diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 3df7fe7..6cf308b 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -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