From bf77863d6cd25757892bebe13367e5d00c6d1bb9 Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 11 Sep 2024 20:35:07 +0200 Subject: [PATCH] QP search for SRG done --- src/GW/GG0W0.f90 | 2 +- src/GW/GGW_QP_graph.f90 | 21 ++++++-- src/GW/{GGW_SigC.f90 => GGW_Re_SigC.f90} | 8 +-- src/GW/{UGW_dSigC.f90 => GGW_Re_dSigC.f90} | 8 +-- ...RGW_SRG_ReSigC.f90 => GGW_SRG_Re_SigC.f90} | 12 ++--- src/GW/GGW_SRG_Re_dSigC.f90 | 53 +++++++++++++++++++ src/GW/RG0W0.f90 | 2 +- src/GW/{RGW_ImSigC.f90 => RGW_Im_SigC.f90} | 0 src/GW/{RGW_ImdSigC.f90 => RGW_Im_dSigC.f90} | 0 src/GW/RGW_QP_graph.f90 | 21 ++++++-- src/GW/{RGW_ReSigC.f90 => RGW_Re_SigC.f90} | 8 +-- src/GW/{RGW_RedSigC.f90 => RGW_Re_dSigC.f90} | 8 +-- src/GW/RGW_SRG_Re_SigC.f90 | 52 ++++++++++++++++++ src/GW/RGW_SRG_Re_dSigC.f90 | 53 +++++++++++++++++++ src/GW/RGW_SRG_self_energy_diag.f90 | 2 +- src/GW/UG0W0.f90 | 2 +- src/GW/UGW_QP_graph.f90 | 21 ++++++-- src/GW/{UGW_SigC.f90 => UGW_Re_SigC.f90} | 8 +-- src/GW/{GGW_dSigC.f90 => UGW_Re_dSigC.f90} | 8 +-- src/GW/UGW_SRG_Re_SigC.f90 | 52 ++++++++++++++++++ src/GW/UGW_SRG_Re_dSigC.f90 | 53 +++++++++++++++++++ src/GW/evGGW.f90 | 2 +- src/GW/evRGW.f90 | 2 +- src/GW/evUGW.f90 | 2 +- 24 files changed, 351 insertions(+), 49 deletions(-) rename src/GW/{GGW_SigC.f90 => GGW_Re_SigC.f90} (82%) rename src/GW/{UGW_dSigC.f90 => GGW_Re_dSigC.f90} (80%) rename src/GW/{RGW_SRG_ReSigC.f90 => GGW_SRG_Re_SigC.f90} (75%) create mode 100644 src/GW/GGW_SRG_Re_dSigC.f90 rename src/GW/{RGW_ImSigC.f90 => RGW_Im_SigC.f90} (100%) rename src/GW/{RGW_ImdSigC.f90 => RGW_Im_dSigC.f90} (100%) rename src/GW/{RGW_ReSigC.f90 => RGW_Re_SigC.f90} (83%) rename src/GW/{RGW_RedSigC.f90 => RGW_Re_dSigC.f90} (80%) create mode 100644 src/GW/RGW_SRG_Re_SigC.f90 create mode 100644 src/GW/RGW_SRG_Re_dSigC.f90 rename src/GW/{UGW_SigC.f90 => UGW_Re_SigC.f90} (82%) rename src/GW/{GGW_dSigC.f90 => UGW_Re_dSigC.f90} (80%) create mode 100644 src/GW/UGW_SRG_Re_SigC.f90 create mode 100644 src/GW/UGW_SRG_Re_dSigC.f90 diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 109b2b0..fea0d23 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -140,7 +140,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) + call GGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) end if diff --git a/src/GW/GGW_QP_graph.f90 b/src/GW/GGW_QP_graph.f90 index 1d82545..6c90770 100644 --- a/src/GW/GGW_QP_graph.f90 +++ b/src/GW/GGW_QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) +subroutine GGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) ! Compute the graphical solution of the QP equation @@ -14,7 +14,9 @@ subroutine GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) integer,intent(in) :: nR integer,intent(in) :: nS + logical,intent(in) :: doSRG double precision,intent(in) :: eta + double precision,intent(in) :: flow double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) @@ -28,7 +30,8 @@ subroutine GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) integer :: nIt integer,parameter :: maxIt = 64 double precision,parameter :: thresh = 1d-6 - double precision,external :: GGW_SigC,GGW_dSigC + double precision,external :: GGW_Re_SigC,GGW_Re_dSigC + double precision,external :: GGW_SRG_Re_SigC,GGW_SRG_Re_dSigC double precision :: SigC,dSigC double precision :: f,df double precision :: w @@ -54,8 +57,18 @@ subroutine GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) nIt = nIt + 1 - SigC = GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) - dSigC = GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + if(doSRG) then + + SigC = GGW_SRG_Re_SigC(p,w,flow,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = GGW_SRG_Re_dSigC(p,w,flow,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + + else + + SigC = GGW_Re_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = GGW_Re_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + + end if + f = w - eHF(p) - SigC df = 1d0/(1d0 - dSigC) diff --git a/src/GW/GGW_SigC.f90 b/src/GW/GGW_Re_SigC.f90 similarity index 82% rename from src/GW/GGW_SigC.f90 rename to src/GW/GGW_Re_SigC.f90 index 6372537..4b9ac88 100644 --- a/src/GW/GGW_SigC.f90 +++ b/src/GW/GGW_Re_SigC.f90 @@ -1,4 +1,4 @@ -double precision function GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function GGW_Re_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute diagonal of the correlation part of the self-energy @@ -27,7 +27,7 @@ double precision function GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - GGW_SigC = 0d0 + GGW_Re_SigC = 0d0 ! Occupied part of the correlation self-energy @@ -35,7 +35,7 @@ double precision function GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(i) + Om(m) num = rho(p,i,m)**2 - GGW_SigC = GGW_SigC + num*eps/(eps**2 + eta**2) + GGW_Re_SigC = GGW_Re_SigC + num*eps/(eps**2 + eta**2) end do end do @@ -45,7 +45,7 @@ double precision function GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(a) - Om(m) num = rho(p,a,m)**2 - GGW_SigC = GGW_SigC + num*eps/(eps**2 + eta**2) + GGW_Re_SigC = GGW_Re_SigC + num*eps/(eps**2 + eta**2) end do end do diff --git a/src/GW/UGW_dSigC.f90 b/src/GW/GGW_Re_dSigC.f90 similarity index 80% rename from src/GW/UGW_dSigC.f90 rename to src/GW/GGW_Re_dSigC.f90 index 8a82fc7..258eb57 100644 --- a/src/GW/UGW_dSigC.f90 +++ b/src/GW/GGW_Re_dSigC.f90 @@ -1,4 +1,4 @@ -double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function GGW_Re_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute the derivative of the correlation part of the self-energy @@ -27,7 +27,7 @@ double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - UGW_dSigC = 0d0 + GGW_Re_dSigC = 0d0 ! Occupied part of the correlation self-energy @@ -35,7 +35,7 @@ double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(i) + Om(m) num = rho(p,i,m)**2 - UGW_dSigC = UGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + GGW_Re_dSigC = GGW_Re_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -45,7 +45,7 @@ double precision function UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(a) - Om(m) num = rho(p,a,m)**2 - UGW_dSigC = UGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + GGW_Re_dSigC = GGW_Re_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do diff --git a/src/GW/RGW_SRG_ReSigC.f90 b/src/GW/GGW_SRG_Re_SigC.f90 similarity index 75% rename from src/GW/RGW_SRG_ReSigC.f90 rename to src/GW/GGW_SRG_Re_SigC.f90 index e15f607..7e94765 100644 --- a/src/GW/RGW_SRG_ReSigC.f90 +++ b/src/GW/GGW_SRG_Re_SigC.f90 @@ -1,4 +1,4 @@ -double precision function RGW_SRG_ReSigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function GGW_SRG_Re_SigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute diagonal of the correlation part of the self-energy @@ -27,15 +27,15 @@ double precision function RGW_SRG_ReSigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - RGW_SRG_ReSigC = 0d0 + GGW_SRG_Re_SigC = 0d0 ! Occupied part of the correlation self-energy do i=nC+1,nO do m=1,nS Dpim = w - e(i) + Om(m) - RGW_SRG_ReSigC = RGW_SRG_ReSigC & - + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim + GGW_SRG_Re_SigC = GGW_SRG_Re_SigC & + + rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim end do end do @@ -44,8 +44,8 @@ double precision function RGW_SRG_ReSigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) do a=nO+1,nBas-nR do m=1,nS Dpam = w - e(a) - Om(m) - RGW_SRG_ReSigC = RGW_SRG_ReSigC & - + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam + GGW_SRG_Re_SigC = GGW_SRG_Re_SigC & + + rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam end do end do diff --git a/src/GW/GGW_SRG_Re_dSigC.f90 b/src/GW/GGW_SRG_Re_dSigC.f90 new file mode 100644 index 0000000..bc04469 --- /dev/null +++ b/src/GW/GGW_SRG_Re_dSigC.f90 @@ -0,0 +1,53 @@ +double precision function GGW_SRG_Re_dSigC(p,w,s,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) + +! Compute the derivative of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: s + double precision,intent(in) :: eta + 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) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: Dpim,Dpam + +! Initialize + + GGW_SRG_Re_dSigC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do m=1,nS + Dpim = w - e(i) + Om(m) + GGW_SRG_Re_dSigC = GGW_SRG_Re_dSigC & + - rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2 + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do m=1,nS + Dpam = w - e(a) - Om(m) + GGW_SRG_Re_dSigC = GGW_SRG_Re_dSigC & + - rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + +end function diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index 32c2a61..272c2d9 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -146,7 +146,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call RGW_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) + call RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) end if diff --git a/src/GW/RGW_ImSigC.f90 b/src/GW/RGW_Im_SigC.f90 similarity index 100% rename from src/GW/RGW_ImSigC.f90 rename to src/GW/RGW_Im_SigC.f90 diff --git a/src/GW/RGW_ImdSigC.f90 b/src/GW/RGW_Im_dSigC.f90 similarity index 100% rename from src/GW/RGW_ImdSigC.f90 rename to src/GW/RGW_Im_dSigC.f90 diff --git a/src/GW/RGW_QP_graph.f90 b/src/GW/RGW_QP_graph.f90 index d7d1ae0..0bbca84 100644 --- a/src/GW/RGW_QP_graph.f90 +++ b/src/GW/RGW_QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) +subroutine RGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) ! Compute the graphical solution of the QP equation @@ -14,7 +14,9 @@ subroutine RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) integer,intent(in) :: nR integer,intent(in) :: nS + logical,intent(in) :: doSRG double precision,intent(in) :: eta + double precision,intent(in) :: flow double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) @@ -28,7 +30,8 @@ subroutine RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) integer :: nIt integer,parameter :: maxIt = 64 double precision,parameter :: thresh = 1d-6 - double precision,external :: RGW_ReSigC,RGW_RedSigC + double precision,external :: RGW_Re_SigC,RGW_Re_dSigC + double precision,external :: RGW_SRG_Re_SigC,RGW_SRG_Re_dSigC double precision :: SigC,dSigC double precision :: f,df double precision :: w @@ -54,8 +57,18 @@ subroutine RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) nIt = nIt + 1 - SigC = RGW_ReSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) - dSigC = RGW_RedSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + if(doSRG) then + + SigC = RGW_SRG_Re_SigC(p,w,flow,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = RGW_SRG_Re_dSigC(p,w,flow,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + + else + + SigC = RGW_Re_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = RGW_Re_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + + end if + f = w - eHF(p) - SigC df = 1d0/(1d0 - dSigC) w = w - df*f diff --git a/src/GW/RGW_ReSigC.f90 b/src/GW/RGW_Re_SigC.f90 similarity index 83% rename from src/GW/RGW_ReSigC.f90 rename to src/GW/RGW_Re_SigC.f90 index 0a565aa..1711027 100644 --- a/src/GW/RGW_ReSigC.f90 +++ b/src/GW/RGW_Re_SigC.f90 @@ -1,4 +1,4 @@ -double precision function RGW_ReSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function RGW_Re_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute diagonal of the correlation part of the self-energy @@ -27,7 +27,7 @@ double precision function RGW_ReSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - RGW_ReSigC = 0d0 + RGW_Re_SigC = 0d0 ! Occupied part of the correlation self-energy @@ -35,7 +35,7 @@ double precision function RGW_ReSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(i) + Om(m) num = 2d0*rho(p,i,m)**2 - RGW_ReSigC = RGW_ReSigC + num*eps/(eps**2 + eta**2) + RGW_Re_SigC = RGW_Re_SigC + num*eps/(eps**2 + eta**2) end do end do @@ -45,7 +45,7 @@ double precision function RGW_ReSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(a) - Om(m) num = 2d0*rho(p,a,m)**2 - RGW_ReSigC = RGW_ReSigC + num*eps/(eps**2 + eta**2) + RGW_Re_SigC = RGW_Re_SigC + num*eps/(eps**2 + eta**2) end do end do diff --git a/src/GW/RGW_RedSigC.f90 b/src/GW/RGW_Re_dSigC.f90 similarity index 80% rename from src/GW/RGW_RedSigC.f90 rename to src/GW/RGW_Re_dSigC.f90 index 9c59415..cabeb6a 100644 --- a/src/GW/RGW_RedSigC.f90 +++ b/src/GW/RGW_Re_dSigC.f90 @@ -1,4 +1,4 @@ -double precision function RGW_RedSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function RGW_Re_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute the derivative of the correlation part of the self-energy @@ -27,7 +27,7 @@ double precision function RGW_RedSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - RGW_RedSigC = 0d0 + RGW_Re_dSigC = 0d0 ! Occupied part of the correlation self-energy @@ -35,7 +35,7 @@ double precision function RGW_RedSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(i) + Om(m) num = 2d0*rho(p,i,m)**2 - RGW_RedSigC = RGW_RedSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + RGW_Re_dSigC = RGW_Re_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -45,7 +45,7 @@ double precision function RGW_RedSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(a) - Om(m) num = 2d0*rho(p,a,m)**2 - RGW_RedSigC = RGW_RedSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + RGW_Re_dSigC = RGW_Re_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do diff --git a/src/GW/RGW_SRG_Re_SigC.f90 b/src/GW/RGW_SRG_Re_SigC.f90 new file mode 100644 index 0000000..f654139 --- /dev/null +++ b/src/GW/RGW_SRG_Re_SigC.f90 @@ -0,0 +1,52 @@ +double precision function RGW_SRG_Re_SigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: s + 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) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: Dpim,Dpam + +! Initialize + + RGW_SRG_Re_SigC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do m=1,nS + Dpim = w - e(i) + Om(m) + RGW_SRG_Re_SigC = RGW_SRG_Re_SigC & + + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do m=1,nS + Dpam = w - e(a) - Om(m) + RGW_SRG_Re_SigC = RGW_SRG_Re_SigC & + + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam + end do + end do + +end function diff --git a/src/GW/RGW_SRG_Re_dSigC.f90 b/src/GW/RGW_SRG_Re_dSigC.f90 new file mode 100644 index 0000000..7d32a3e --- /dev/null +++ b/src/GW/RGW_SRG_Re_dSigC.f90 @@ -0,0 +1,53 @@ +double precision function RGW_SRG_Re_dSigC(p,w,s,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) + +! Compute the derivative of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: s + double precision,intent(in) :: eta + 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) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: Dpim,Dpam + +! Initialize + + RGW_SRG_Re_dSigC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do m=1,nS + Dpim = w - e(i) + Om(m) + RGW_SRG_Re_dSigC = RGW_SRG_Re_dSigC & + - 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2 + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do m=1,nS + Dpam = w - e(a) - Om(m) + RGW_SRG_Re_dSigC = RGW_SRG_Re_dSigC & + - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + +end function diff --git a/src/GW/RGW_SRG_self_energy_diag.f90 b/src/GW/RGW_SRG_self_energy_diag.f90 index eef9f29..4a78b74 100644 --- a/src/GW/RGW_SRG_self_energy_diag.f90 +++ b/src/GW/RGW_SRG_self_energy_diag.f90 @@ -110,7 +110,7 @@ subroutine RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM, do a=nO+1,nOrb-nR do m=1,nS Dpam = e(p) - e(a) - Om(m) - Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 + Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 end do end do end do diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 0975f97..9a22243 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -164,7 +164,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD if(is==1) write(*,*)' Spin-up orbitals ' if(is==2) write(*,*)' Spin-down orbitals ' - call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is), & + call UGW_QP_graph(doSRG,eta,flow,nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is), & Om,rho(:,:,:,is),eGWlin(:,is),eHF(:,is),eGW(:,is),Z(:,is)) end do diff --git a/src/GW/UGW_QP_graph.f90 b/src/GW/UGW_QP_graph.f90 index f9606d5..1da8b86 100644 --- a/src/GW/UGW_QP_graph.f90 +++ b/src/GW/UGW_QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) +subroutine UGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) ! Compute the graphical solution of the QP equation @@ -14,7 +14,9 @@ subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) integer,intent(in) :: nR integer,intent(in) :: nS + logical,intent(in) :: doSRG double precision,intent(in) :: eta + double precision,intent(in) :: flow double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS,nspin) @@ -28,7 +30,8 @@ subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) integer :: nIt integer,parameter :: maxIt = 64 double precision,parameter :: thresh = 1d-6 - double precision,external :: UGW_SigC,UGW_dSigC + double precision,external :: UGW_Re_SigC,UGW_Re_dSigC + double precision,external :: UGW_SRG_Re_SigC,UGW_SRG_Re_dSigC double precision :: SigC,dSigC double precision :: f,df double precision :: w @@ -54,8 +57,18 @@ subroutine UGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) nIt = nIt + 1 - SigC = UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) - dSigC = UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + if(doSRG) then + + SigC = UGW_SRG_Re_SigC(p,w,flow,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = UGW_SRG_Re_dSigC(p,w,flow,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + + else + + SigC = UGW_Re_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + dSigC = UGW_Re_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) + + end if + f = w - eHF(p) - SigC df = 1d0/(1d0 - dSigC) diff --git a/src/GW/UGW_SigC.f90 b/src/GW/UGW_Re_SigC.f90 similarity index 82% rename from src/GW/UGW_SigC.f90 rename to src/GW/UGW_Re_SigC.f90 index 488c983..b708884 100644 --- a/src/GW/UGW_SigC.f90 +++ b/src/GW/UGW_Re_SigC.f90 @@ -1,4 +1,4 @@ -double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function UGW_Re_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute diagonal of the correlation part of the self-energy @@ -27,7 +27,7 @@ double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - UGW_SigC = 0d0 + UGW_Re_SigC = 0d0 ! Occupied part of the correlation self-energy @@ -35,7 +35,7 @@ double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(i) + Om(m) num = rho(p,i,m)**2 - UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2) + UGW_Re_SigC = UGW_Re_SigC + num*eps/(eps**2 + eta**2) end do end do @@ -45,7 +45,7 @@ double precision function UGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(a) - Om(m) num = rho(p,a,m)**2 - UGW_SigC = UGW_SigC + num*eps/(eps**2 + eta**2) + UGW_Re_SigC = UGW_Re_SigC + num*eps/(eps**2 + eta**2) end do end do diff --git a/src/GW/GGW_dSigC.f90 b/src/GW/UGW_Re_dSigC.f90 similarity index 80% rename from src/GW/GGW_dSigC.f90 rename to src/GW/UGW_Re_dSigC.f90 index e1bd2c4..cb113bf 100644 --- a/src/GW/GGW_dSigC.f90 +++ b/src/GW/UGW_Re_dSigC.f90 @@ -1,4 +1,4 @@ -double precision function GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) +double precision function UGW_Re_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Compute the derivative of the correlation part of the self-energy @@ -27,7 +27,7 @@ double precision function GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - GGW_dSigC = 0d0 + UGW_Re_dSigC = 0d0 ! Occupied part of the correlation self-energy @@ -35,7 +35,7 @@ double precision function GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(i) + Om(m) num = rho(p,i,m)**2 - GGW_dSigC = GGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + UGW_Re_dSigC = UGW_Re_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -45,7 +45,7 @@ double precision function GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) do m=1,nS eps = w - e(a) - Om(m) num = rho(p,a,m)**2 - GGW_dSigC = GGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + UGW_Re_dSigC = UGW_Re_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do diff --git a/src/GW/UGW_SRG_Re_SigC.f90 b/src/GW/UGW_SRG_Re_SigC.f90 new file mode 100644 index 0000000..e9c8d3c --- /dev/null +++ b/src/GW/UGW_SRG_Re_SigC.f90 @@ -0,0 +1,52 @@ +double precision function UGW_SRG_Re_SigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: s + 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) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: Dpim,Dpam + +! Initialize + + UGW_SRG_Re_SigC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do m=1,nS + Dpim = w - e(i) + Om(m) + UGW_SRG_Re_SigC = UGW_SRG_Re_SigC & + + rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do m=1,nS + Dpam = w - e(a) - Om(m) + UGW_SRG_Re_SigC = UGW_SRG_Re_SigC & + + rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam + end do + end do + +end function diff --git a/src/GW/UGW_SRG_Re_dSigC.f90 b/src/GW/UGW_SRG_Re_dSigC.f90 new file mode 100644 index 0000000..0590dc4 --- /dev/null +++ b/src/GW/UGW_SRG_Re_dSigC.f90 @@ -0,0 +1,53 @@ +double precision function UGW_SRG_Re_dSigC(p,w,s,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) + +! Compute the derivative of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: s + double precision,intent(in) :: eta + 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) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: Dpim,Dpam + +! Initialize + + UGW_SRG_Re_dSigC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do m=1,nS + Dpim = w - e(i) + Om(m) + UGW_SRG_Re_dSigC = UGW_SRG_Re_dSigC & + - rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2 + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do m=1,nS + Dpam = w - e(a) - Om(m) + UGW_SRG_Re_dSigC = UGW_SRG_Re_dSigC & + - rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + +end function diff --git a/src/GW/evGGW.f90 b/src/GW/evGGW.f90 index f8c59f0..1ef3bf6 100644 --- a/src/GW/evGGW.f90 +++ b/src/GW/evGGW.f90 @@ -146,7 +146,7 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) + call GGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) end if diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index 438ffb4..9455b23 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -151,7 +151,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call RGW_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) + call RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) end if diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 80f4ce0..4d7d44a 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -181,7 +181,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE if(is==1) write(*,*)' Spin-up orbitals ' if(is==2) write(*,*)' Spin-down orbitals ' - call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is), & + call UGW_QP_graph(doSRG,eta,flow,nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is), & Om,rho(:,:,:,is),eOld(:,is),eOld(:,is),eGW(:,is),Z(:,is)) end do