SRG for MOM U self energy

Author: lburth

src/GW/CVS_UGW_SRG_self_energy_diag.f90

@@ -1,4 +1,4 @@
-subroutine CVS_UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nSt,nCVS,nFC,occupations,virtuals,e,Om,rho,EcGM,Sig,Z)
+subroutine CVS_UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nSt,nCVS,nFC,occupations,virtuals,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -24,24 +24,65 @@ subroutine CVS_UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nSt,nCVS,nFC,occup
 
 ! Local variables
 
-  integer                       :: i,a,p,m
+  integer                       :: i,a,p,m,ispin
   double precision              :: num,eps
   double precision              :: s
   double precision              :: Dpim,Dpam,Diam
 
 ! Output variables
 
-  double precision,intent(out)  :: Sig(nBas,nspin)
+  double precision,intent(out)  :: SigC(nBas,nspin)
   double precision,intent(out)  :: Z(nBas,nspin)
   double precision              :: EcGM(nspin)
 
 ! Initialize 
 
-  Sig(:,:) = 0d0
+  SigC(:,:) = 0d0
   Z(:,:)   = 0d0
   EcGM(:)  = 0d0
   s = flow
 
-  print *, "SRG not implemented for CVS/MOM yet"
 
+  do ispin=1,nspin
+    ! Occupied part of the correlation self-energy
+
+    do p=1,nBas
+      do i=1,nO(ispin)-nFC(ispin)
+        do m=1,nSt
+           Dpim = e(p,ispin) - e(occupations(i,ispin),ispin) + Om(m)
+           SigC(p,ispin) = SigC(p,ispin) + rho(p,occupations(i,ispin),m,ispin)**2&
+                                         * (1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
+        end do
+      end do
+    end do
+
+    ! Virtual part of the correlation self-energy
+
+    do p=1,nBas
+      do a=1+nCVS(ispin),nBas - nO(ispin)
+        do m=1,nSt
+          Dpam = e(p,ispin) - e(virtuals(a,ispin),ispin) - Om(m)
+          SigC(p,ispin) = SigC(p,ispin) + rho(p,virtuals(a,ispin),m,ispin)**2&
+                                        * (1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
+        end do
+      end do
+    end do
+
+    ! GM correlation energy
+
+    do i=1,nO(ispin)-nFC(ispin)
+      do a=nCVS(ispin)+1,nBas-nO(ispin)
+        do m=1,nSt
+          Diam = e(virtuals(a,ispin),ispin) - e(occupations(i,ispin),ispin) + Om(m)
+          EcGM = EcGM - rho(virtuals(a,ispin),occupations(i,ispin),m,ispin)&
+                      * rho(virtuals(a,ispin),occupations(i,ispin),m,ispin)&
+                      * (1d0-exp(-2d0*s*Diam*Diam))/Diam 
+        end do
+      end do
+    end do
+  end do
+
+! Compute renormalization factor from derivative 
+
+  Z(:,:) = 1d0/(1d0 - Z(:,:))
 end subroutine