From e19e8708053d0b2e9af04a04359e00c09a109746 Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 24 Aug 2023 11:46:45 +0200 Subject: [PATCH] fix root search in GTpp and plot routine --- input/methods | 4 +- src/GT/G0T0pp.f90 | 3 + src/GT/GTeh_QP_graph.f90 | 3 +- src/GT/GTpp_QP_graph.f90 | 24 +++++--- src/GT/GTpp_plot_self_energy.f90 | 102 +++++++++++++++++++++++++++++++ 5 files changed, 123 insertions(+), 13 deletions(-) create mode 100644 src/GT/GTpp_plot_self_energy.f90 diff --git a/input/methods b/input/methods index ddf1e5d..91124e8 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW - T F F F F F + F F F F F F # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh - F F F T F F + T F F F F F # * unrestricted version available diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 5a7384e..63c0ddb 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -197,6 +197,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp end if + call GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t) + !---------------------------------------------- ! Dump results !---------------------------------------------- diff --git a/src/GT/GTeh_QP_graph.f90 b/src/GT/GTeh_QP_graph.f90 index 4e22971..4db7b04 100644 --- a/src/GT/GTeh_QP_graph.f90 +++ b/src/GT/GTeh_QP_graph.f90 @@ -1,5 +1,6 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) + ! Compute the graphical solution of the QP equation implicit none @@ -67,8 +68,8 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) if(nIt == maxIt) then - write(*,*) 'Newton root search has not converged!' eGT(p) = eGTlin(p) + write(*,*) 'Newton root search has not converged!' else diff --git a/src/GT/GTpp_QP_graph.f90 b/src/GT/GTpp_QP_graph.f90 index 13aa90b..a5dc219 100644 --- a/src/GT/GTpp_QP_graph.f90 +++ b/src/GT/GTpp_QP_graph.f90 @@ -1,10 +1,13 @@ 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,Z) +! Compute the graphical solution of the QP equation + implicit none include 'parameters.h' -! Iput variables +! Input variables + integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -23,6 +26,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s double precision,intent(in) :: eGTlin(nBas) ! Local variables + integer :: p integer :: nIt integer,parameter :: maxIt = 64 @@ -41,6 +45,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s dSigC = 0d0 ! Run Newton's algorithm to find the root + do p=nC+1,nBas-nR write(*,*) '-----------------' @@ -48,23 +53,21 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s write(*,*) '-----------------' w = eGTlin(p) - write(*,*) 'HERE', eGTlin(p), eHF(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 + 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) - f = w - eHF(p) - SigC - df = 1d0/(1d0 - dSigC) - - w = w - f/df + SigC = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGTlin,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t) + dSigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGTlin,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t) + 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 + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f end do @@ -80,6 +83,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV write(*,*) + end if diff --git a/src/GT/GTpp_plot_self_energy.f90 b/src/GT/GTpp_plot_self_energy.f90 new file mode 100644 index 0000000..62f418d --- /dev/null +++ b/src/GT/GTpp_plot_self_energy.f90 @@ -0,0 +1,102 @@ +subroutine GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t) + +! Dump several GTpp quantities for external plotting + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt + + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: Om1s(nVVs),Om1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2s(nOOs),Om2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) + +! Local variables + + integer :: p,g + integer :: nGrid + double precision :: wmin,wmax,dw + double precision,external :: GTpp_SigC,GTpp_dSigC + double precision,allocatable :: w(:) + double precision,allocatable :: SigC(:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: S(:,:) + +! Construct grid + + nGrid = 1000 + allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid)) + +! Initialize + + SigC(:,:) = 0d0 + Z(:,:) = 0d0 + +! Minimum and maximum frequency values + + wmin = -5d0 + wmax = +5d0 + dw = (wmax - wmin)/dble(ngrid) + + do g=1,nGrid + w(g) = wmin + dble(g)*dw + enddo + +! Occupied part of the self-energy and renormalization factor + + do g=1,nGrid + do p=nC+1,nBas-nR + + SigC(p,g) = GTpp_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t) + Z(p,g) = GTpp_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t) + + end do + end do + + Z(:,:) = 1d0/(1d0 + Z(:,:)) + +! Compute spectral function + + do g=1,nGrid + do p=nC+1,nBas-nR + S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2) + enddo + enddo + + S(:,:) = S(:,:)/pi + +! Dump quantities in files as a function of w + + open(unit=8 ,file='GTpp_SigC.dat') + open(unit=9 ,file='GTpp_freq.dat') + open(unit=10 ,file='GTpp_Z.dat') + open(unit=11 ,file='GTpp_A.dat') + + do g=1,nGrid + write(8 ,*) w(g)*HaToeV,(SigC(p,g)*HaToeV,p=nC+1,nBas-nR) + write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR) + write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR) + write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR) + enddo + +! Closing files + + close(unit=8) + close(unit=9) + close(unit=10) + close(unit=11) + +end subroutine