diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90
index a5b3c0f..63b5453 100644
--- a/src/GF/RG0F2.f90
+++ b/src/GF/RG0F2.f90
@@ -71,13 +71,15 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
 
   if(doSRG) then 
 
-    call RGF2_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
+    call RGF2_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,eHF,ERI,Ec,SigC,Z)
 
   else
 
-    call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
+    call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,Ec,SigC,Z)
 
   end if
+
+  print*,Ec
   
   eGFlin(:) = eHF(:) + Z(:)*SigC(:)
 
@@ -98,7 +100,6 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
 
   ! Print results
 
-  call RMP2(.false.,doSRG,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
   call print_RG0F2(nOrb,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
 
 ! Perform BSE@GF2 calculation
diff --git a/src/GF/RGF2_SRG_self_energy.f90 b/src/GF/RGF2_SRG_self_energy.f90
index fd48d95..9aed51f 100644
--- a/src/GF/RGF2_SRG_self_energy.f90
+++ b/src/GF/RGF2_SRG_self_energy.f90
@@ -1,4 +1,4 @@
-subroutine RGF2_SRG_self_energy(flow,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_SRG_self_energy(flow,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute GF2 self-energy and its renormalization factor
 
@@ -21,13 +21,14 @@ subroutine RGF2_SRG_self_energy(flow,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
   integer                       :: i,j,a,b
   integer                       :: p,q
   double precision              :: eps_p,eps_q
-  double precision              :: num
+  double precision              :: num,eps
 
   double precision              :: s
   double precision              :: kappa
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas,nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -121,4 +122,26 @@ subroutine RGF2_SRG_self_energy(flow,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
   Z(:) = 1d0/(1d0 - Z(:))
 
+!----------------------------!
+! Compute correlation energy !
+!----------------------------!
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          kappa = 1d0 - exp(-2d0*s*eps**2)
+          num = kappa*(2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+
 end subroutine
diff --git a/src/GF/RGF2_SRG_self_energy_diag.f90 b/src/GF/RGF2_SRG_self_energy_diag.f90
index 2360b42..d58f5a2 100644
--- a/src/GF/RGF2_SRG_self_energy_diag.f90
+++ b/src/GF/RGF2_SRG_self_energy_diag.f90
@@ -1,4 +1,4 @@
-subroutine RGF2_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute diagonal part of the GF2 self-energy and its renormalization factor
 
@@ -28,6 +28,7 @@ subroutine RGF2_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -82,4 +83,26 @@ subroutine RGF2_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
   Z(:) = 1d0/(1d0 - Z(:))
 
+!----------------------------!  
+! Compute correlation energy !  
+!----------------------------!  
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          kappa = 1d0 - exp(-2d0*s*eps**2)
+          num = kappa*(2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+
 end subroutine 
diff --git a/src/GF/RGF2_self_energy.f90 b/src/GF/RGF2_self_energy.f90
index ec09c0b..f4d3a34 100644
--- a/src/GF/RGF2_self_energy.f90
+++ b/src/GF/RGF2_self_energy.f90
@@ -1,4 +1,4 @@
-subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute GF2 self-energy and its renormalization factor
 
@@ -25,6 +25,7 @@ subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas,nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -73,4 +74,23 @@ subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
   Z(:) = 1d0/(1d0 - Z(:))
 
+! Compute MP2 correlation energy
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          num = (2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+  
 end subroutine 
diff --git a/src/GF/RGF2_self_energy_diag.f90 b/src/GF/RGF2_self_energy_diag.f90
index 71822ab..977dcf3 100644
--- a/src/GF/RGF2_self_energy_diag.f90
+++ b/src/GF/RGF2_self_energy_diag.f90
@@ -1,4 +1,4 @@
-subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute diagonal part of the GF2 self-energy and its renormalization factor
 
@@ -25,6 +25,7 @@ subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -66,6 +67,26 @@ subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
       end do
     end do
   end do
+
   Z(:) = 1d0/(1d0 - Z(:))
 
+! Compute correlation energy
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          num = (2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+
 end subroutine 
diff --git a/src/GF/evRGF2.f90 b/src/GF/evRGF2.f90
index 2a38262..60d8d7e 100644
--- a/src/GF/evRGF2.f90
+++ b/src/GF/evRGF2.f90
@@ -96,11 +96,11 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
     if(doSRG) then 
 
-      call RGF2_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z)
+      call RGF2_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,eGF,ERI,Ec,SigC,Z)
 
     else
 
-      call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z)
+      call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,Ec,SigC,Z)
 
     end if
 
@@ -121,7 +121,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
     ! Print results
 
-    call RMP2(.false.,doSRG,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
     call print_evRGF2(nOrb,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
 
     ! DIIS extrapolation
diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90
index ff261b3..499a3e9 100644
--- a/src/GF/qsRGF2.f90
+++ b/src/GF/qsRGF2.f90
@@ -177,14 +177,16 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
 
     if(doSRG) then
 
-      call RGF2_SRG_self_energy(flow,nOrb,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z)
+      call RGF2_SRG_self_energy(flow,nOrb,nC,nO,nV,nR,eGF,ERI_MO,Ec,SigC,Z)
 
     else
 
-      call RGF2_self_energy(eta,nOrb,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z)
+      call RGF2_self_energy(eta,nOrb,nC,nO,nV,nR,eGF,ERI_MO,Ec,SigC,Z)
 
     end if
 
+    print*,Ec
+
     ! Make correlation self-energy Hermitian and transform it back to AO basis
    
     SigC = 0.5d0*(SigC + transpose(SigC))
@@ -221,10 +223,6 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
 
     Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
 
-    ! Correlation energy
-
-    call RMP2(.false.,doSRG,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec)
-
     ! Total energy
 
     EqsGF2 = ET + EV + EJ + Ex + Ec