QP search for SRG done

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
 

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)
+      if(doSRG) then
-      dSigC = GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
+
+        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)
     

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
 

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
 

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 &
+      GGW_SRG_Re_SigC = GGW_SRG_Re_SigC &
-                     + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
+                      + 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 & 
+      GGW_SRG_Re_SigC = GGW_SRG_Re_SigC & 
-                     + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
+                      + rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
     end do
   end do
 

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 

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
 

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)
+      if(doSRG) then
-      dSigC = RGW_RedSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
+
+        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

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
 

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
 

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 

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 

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

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
  

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)
+      if(doSRG) then
-      dSigC = UGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
+
+        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)
     

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
 

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
 

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 

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 

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
 

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
 

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