From 990cf8462bcd7494da41b33a83005e8317742884 Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 30 Jan 2024 12:06:21 +0100 Subject: [PATCH] cleaning up GW+C --- src/GW/GW_plot_self_energy.f90 | 8 +--- src/GW/RG0W0.f90 | 4 +- src/GW/RGWC.f90 | 83 ++++++++++++++++++---------------- 3 files changed, 49 insertions(+), 46 deletions(-) diff --git a/src/GW/GW_plot_self_energy.f90 b/src/GW/GW_plot_self_energy.f90 index d668700..85b36da 100644 --- a/src/GW/GW_plot_self_energy.f90 +++ b/src/GW/GW_plot_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) +subroutine GW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) ! Dump several GW quantities for external plotting @@ -7,6 +7,7 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) ! Input variables + double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -20,7 +21,6 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) ! Local variables - double precision :: eta integer :: p,g integer :: nGrid double precision :: wmin,wmax,dw @@ -30,10 +30,6 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) double precision,allocatable :: Z(:,:) double precision,allocatable :: A(:,:) -! Broadening parameter - - eta = 0.01d0 - ! Construct grid nGrid = 5000 diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index d172855..fbb05d1 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -148,13 +148,13 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Plot self-energy, renormalization factor, and spectral function - call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) + call GW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) !--------------------! ! Cumulant expansion ! !--------------------! - call RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) + call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) ! Compute the RPA correlation energy diff --git a/src/GW/RGWC.f90 b/src/GW/RGWC.f90 index c683438..cfc734c 100644 --- a/src/GW/RGWC.f90 +++ b/src/GW/RGWC.f90 @@ -1,4 +1,4 @@ -subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) +subroutine RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) ! Perform GW+C calculation @@ -9,6 +9,7 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) logical,intent(in) :: dotest + double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -26,17 +27,17 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) integer :: p,q,i,a,m integer :: iSat - double precision,parameter :: cutoff = 0d0 + double precision :: num,eps + double precision,parameter :: cutoff = 1d-2 logical,parameter :: do_hole_branch = .true. logical,parameter :: do_electron_branch = .false. double precision,allocatable :: de(:,:,:) double precision,allocatable :: ZC(:) - double precision,allocatable :: ZSat(:,:) - double precision,allocatable :: eSat(:,:) + double precision,allocatable :: ZSat(:,:,:) + double precision,allocatable :: eSat(:,:,:) - double precision :: eta integer :: g integer :: nGrid double precision :: wmin,wmax,dw @@ -58,7 +59,7 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) ! Memory allocation - allocate(ZC(nBas),eSat(nBas,nS),ZSat(nBas,nS)) + allocate(ZC(nBas),eSat(nBas,nBas,nS),ZSat(nBas,nBas,nS)) ! Useful quantities @@ -67,7 +68,7 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) do p=nC+1,nBas-nR do i=nC+1,nO do m=1,nS - de(p,i,m) = eGW(p) - eGW(i) + Om(m) + de(p,i,m) = eHF(p) - eGW(i) + Om(m) end do end do end do @@ -75,7 +76,7 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) do p=nC+1,nBas-nR do a=nO+1,nBas-nR do m=1,nS - de(p,a,m) = eGW(p) - eGW(a) - Om(m) + de(p,a,m) = eHF(p) - eGW(a) - Om(m) end do end do end do @@ -87,12 +88,14 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) do p=nC+1,nBas-nR do q=nC+1,nBas-nR do m=1,nS - ZC(p) = ZC(p) + 2d0*rho(p,q,m)**2/de(p,q,m)**2 + num = 2d0*rho(p,q,m)**2 + eps = de(p,q,m) + ZC(p) = ZC(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do end do - ZC(:) = exp(-ZC(:)) + ZC(:) = exp(ZC(:)) ! Dump results @@ -115,12 +118,13 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) if(do_hole_branch) then - ZSat(:,:) = 0d0 do i=nC+1,nO - do m=1,nS - eSat(i,m) = eGW(i) - Om(m) - do q=nC+1,nBas-nR - ZSat(i,m) = ZSat(i,m) + ZC(i)*2d0*rho(i,q,m)**2/de(i,q,m)**2 + do q=nC+1,nBas-nR + do m=1,nS + eps = de(i,q,m) + num = ZC(i)*2d0*rho(i,q,m)**2 + eSat(i,q,m) = eGW(i) + eps + ZSat(i,q,m) = num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do end do @@ -133,10 +137,12 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) write(*,*)'-------------------------------------------------------------------------------' iSat = 0 do i=nC+1,nO - do m=1,nS - iSat = iSat + 1 - if(ZSat(i,m) > cutoff) & - write(*,'(1X,I5,1X,I5,1X,I5,F15.6,1X,F15.6,1X)') iSat,i,m,eSat(i,m)*HaToeV,ZSat(i,m) + do q=nC+1,nBas-nR + do m=1,nS + iSat = iSat + 1 + if(ZSat(i,q,m) > cutoff) & + write(*,'(1X,I5,1X,I5,1X,I5,1X,I5,F15.6,1X,F15.6,1X)') iSat,i,q,m,eSat(i,q,m)*HaToeV,ZSat(i,q,m) + end do end do end do write(*,*)'-------------------------------------------------------------------------------' @@ -148,12 +154,13 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) if(do_electron_branch) then - ZSat(:,:) = 0d0 do a=nO+1,nBas-nR - do m=1,nS - eSat(a,m) = eGW(a) + Om(m) - do q=nC+1,nBas-nR - ZSat(a,m) = ZSat(a,m) + ZC(a)*2d0*rho(a,q,m)**2/de(a,q,m)**2 + do q=nC+1,nBas-nR + do m=1,nS + eps = de(a,q,m) + num = ZC(a)*2d0*rho(a,q,m)**2 + eSat(a,q,m) = eGW(a) + eps + ZSat(a,q,m) = num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do end do @@ -166,10 +173,12 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) write(*,*)'-------------------------------------------------------------------------------' iSat = 0 do a=nO+1,nBas-nR - do m=1,nS - iSat = iSat + 1 - if(ZSat(a,m) > cutoff) & - write(*,'(1X,I5,1X,I5,1X,I5,F15.6,1X,F15.6,1X)') iSat,a,m,eSat(a,m)*HaToeV,ZSat(a,m) + do q=nC+1,nBas-nR + do m=1,nS + iSat = iSat + 1 + if(ZSat(a,q,m) > cutoff) & + write(*,'(1X,I5,1X,I5,I5,1X,1X,I5,F15.6,1X,F15.6,1X)') iSat,a,m,eSat(a,q,m)*HaToeV,ZSat(a,q,m) + end do end do end do write(*,*)'-------------------------------------------------------------------------------' @@ -177,10 +186,6 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) end if -! Broadening parameter - - eta = 0.01d0 - ! Construct grid nGrid = 1000 @@ -202,10 +207,8 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) do g=1,nGrid do p=nC+1,nBas-nR - ReSigC(p,g) = GW_ReSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho) ImSigC(p,g) = GW_ImSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho) - end do end do @@ -237,8 +240,10 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) do g=1,nGrid do i=nC+1,nO - do m=1,nS - AGWC(i,g) = AGWC(i,g) + ZSat(i,m)*abs(ImSigC(i,g))/((w(g) - eSat(i,m))**2 + ImSigC(i,g)**2) + do q=nC+1,nBas-nR + do m=1,nS + AGWC(i,g) = AGWC(i,g) + ZSat(i,q,m)*abs(ImSigC(i,g))/((w(g) - eSat(i,q,m))**2 + ImSigC(i,g)**2) + end do end do end do end do @@ -249,8 +254,10 @@ subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z) do g=1,nGrid do a=nO+1,nBas-nR - do m=1,nS - AGWC(a,g) = AGWC(a,g) + ZSat(a,m)*abs(ImSigC(a,g))/((w(g) - eSat(a,m))**2 + ImSigC(a,g)**2) + do q=nC+1,nBas-nR + do m=1,nS + AGWC(a,g) = AGWC(a,g) + ZSat(a,q,m)*abs(ImSigC(a,g))/((w(g) - eSat(a,q,m))**2 + ImSigC(a,g)**2) + end do end do end do end do