From 7efca5605431caf01b3bf0dc831208ac1cdf9d7a Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 24 Aug 2023 00:02:23 +0200 Subject: [PATCH] plotting routines --- input/methods | 4 +- input/options | 2 +- src/GT/G0T0eh.f90 | 2 + src/GT/GTeh_QP_graph.f90 | 25 ++++---- src/GT/GTeh_plot_self_energy.f90 | 98 ++++++++++++++++++++++++++++++++ src/GT/print_G0T0pp.f90 | 6 +- src/GW/G0W0.f90 | 2 + src/GW/GW_QP_graph.f90 | 3 +- src/GW/GW_plot_self_energy.f90 | 97 +++++++++++++++++++++++++++++++ 9 files changed, 221 insertions(+), 18 deletions(-) create mode 100644 src/GT/GTeh_plot_self_energy.f90 create mode 100644 src/GW/GW_plot_self_energy.f90 diff --git a/input/methods b/input/methods index 6c28410..ddf1e5d 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 - F F F F F F + T F F F F F # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh - F F F F F T + F F F T F F # * unrestricted version available diff --git a/input/options b/input/options index 57f23d6..acdfcbb 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg 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.1 F F + 256 0.00001 T 5 F 0.0 F F # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index a876174..5074e0b 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -149,6 +149,8 @@ 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) + ! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) diff --git a/src/GT/GTeh_QP_graph.f90 b/src/GT/GTeh_QP_graph.f90 index bd5dc6e..5a64ae3 100644 --- a/src/GT/GTeh_QP_graph.f90 +++ b/src/GT/GTeh_QP_graph.f90 @@ -1,9 +1,12 @@ 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 include 'parameters.h' -! Iput variables +! Input variables + integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -20,6 +23,7 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) double precision,intent(in) :: eGTlin(nBas) ! Local variables + integer :: p integer :: nIt integer,parameter :: maxIt = 64 @@ -34,10 +38,8 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) double precision,intent(out) :: eGT(nBas) double precision,intent(out) :: Z(nBas) - sigC = 0d0 - dsigC = 0d0 - ! Run Newton's algorithm to find the root + do p=nC+1,nBas-nR write(*,*) '-----------------' @@ -51,16 +53,16 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) do while (abs(f) > thresh .and. nIt < maxIt) - nIt = nIt + 1 + 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 + 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 @@ -73,6 +75,7 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) eGT(p) = w Z(p) = df + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV write(*,*) diff --git a/src/GT/GTeh_plot_self_energy.f90 b/src/GT/GTeh_plot_self_energy.f90 new file mode 100644 index 0000000..f1a07f3 --- /dev/null +++ b/src/GT/GTeh_plot_self_energy.f90 @@ -0,0 +1,98 @@ +subroutine GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR) + +! Dump several GTeh 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) :: nS + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rhoL(nBas,nBas,nS) + double precision,intent(in) :: rhoR(nBas,nBas,nS) + +! Local variables + + integer :: p,g + integer :: nGrid + double precision :: wmin,wmax,dw + double precision,external :: GTeh_SigC,GTeh_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) = GTeh_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) + Z(p,g) = GTeh_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) + + 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='GTeh_SigC.dat') + open(unit=9 ,file='GTeh_freq.dat') + open(unit=10 ,file='GTeh_Z.dat') + open(unit=11 ,file='GTeh_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 diff --git a/src/GT/print_G0T0pp.f90 b/src/GT/print_G0T0pp.f90 index 2cc1ec2..a77a3a5 100644 --- a/src/GT/print_G0T0pp.f90 +++ b/src/GT/print_G0T0pp.f90 @@ -48,9 +48,9 @@ subroutine print_G0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) enddo write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO energy (eV) =',eGT(HOMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp LUMO energy (eV) =',eGT(LUMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO-LUMO gap (eV) =',Gap*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO energy =',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp LUMO energy =',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO-LUMO gap =',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp correlation energy (singlet) =',EcRPA(1),' au' write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp correlation energy (triplet) =',EcRPA(2),' au' diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index cb5917f..9fcadb3 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -140,6 +140,8 @@ 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) + ! Compute the RPA correlation energy call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) diff --git a/src/GW/GW_QP_graph.f90 b/src/GW/GW_QP_graph.f90 index d88b6a9..c37a872 100644 --- a/src/GW/GW_QP_graph.f90 +++ b/src/GW/GW_QP_graph.f90 @@ -13,6 +13,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z) integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS + double precision,intent(in) :: eta double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: Om(nS) @@ -55,7 +56,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z) 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 + f = w - eHF(p) - sigC df = 1d0/(1d0 - dsigC) w = w - df*f diff --git a/src/GW/GW_plot_self_energy.f90 b/src/GW/GW_plot_self_energy.f90 new file mode 100644 index 0000000..0f079f1 --- /dev/null +++ b/src/GW/GW_plot_self_energy.f90 @@ -0,0 +1,97 @@ +subroutine GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) + +! Dump several GW 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) :: nS + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: p,g + integer :: nGrid + double precision :: wmin,wmax,dw + double precision,external :: GW_SigC,GW_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) = GW_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho) + Z(p,g) = GW_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho) + + 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='GW_SigC.dat') + open(unit=9 ,file='GW_freq.dat') + open(unit=10 ,file='GW_Z.dat') + open(unit=11 ,file='GW_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