diff --git a/src/MBPT/G0F2.f90 b/src/MBPT/G0F2.f90
index db22b4b..303e5ea 100644
--- a/src/MBPT/G0F2.f90
+++ b/src/MBPT/G0F2.f90
@@ -30,15 +30,11 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
 
 ! Local variables
 
-  double precision              :: eps
-  double precision              :: V
   double precision              :: EcBSE(nspin)
   double precision,allocatable  :: eGF2(:)
-  double precision,allocatable  :: Sig(:)
+  double precision,allocatable  :: SigC(:)
   double precision,allocatable  :: Z(:)
 
-  integer                       :: i,j,a,b,p
-
 ! Hello world
 
   write(*,*)
@@ -49,7 +45,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
 
 ! Memory allocation
 
-  allocate(Sig(nBas),Z(nBas),eGF2(nBas))
+  allocate(SigC(nBas),Z(nBas),eGF2(nBas))
 
   if(linearize) then 
   
@@ -60,56 +56,21 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
 
 ! Frequency-dependent second-order contribution
 
-  Sig(:) = 0d0
-  Z(:)   = 0d0
-
-  do p=nC+1,nBas-nR
-    do i=nC+1,nO
-      do j=nC+1,nO
-        do a=nO+1,nBas-nR
-
-          eps = eHF(p) + eHF(a) - eHF(i) - eHF(j)
-          V  = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
-
-          Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2)
-          Z(p)   = Z(p)   - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-        end do
-      end do
-    end do
-  end do
-
-  do p=nC+1,nBas-nR
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        do b=nO+1,nBas-nR
-
-          eps = eHF(p) + eHF(i) - eHF(a) - eHF(b)
-          V  = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
-
-          Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2)
-          Z(p)   = Z(p)   - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-        end do
-      end do
-    end do
-  end do
-
-  Z(:) = 1d0/(1d0 - Z(:))
+  call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z)
 
   if(linearize) then
 
-    eGF2(:) = eHF(:) + Z(:)*Sig(:)
+    eGF2(:) = eHF(:) + Z(:)*SigC(:)
 
   else
 
-    eGF2(:) = eHF(:) + Sig(:)
+    eGF2(:) = eHF(:) + SigC(:)
 
   end if
 
   ! Print results
 
-  call print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
+  call print_G0F2(nBas,nO,eHF,SigC,eGF2,Z)
 
 ! Perform BSE2 calculation
 
diff --git a/src/MBPT/evGF2.f90 b/src/MBPT/evGF2.f90
index 14bd86f..ff39619 100644
--- a/src/MBPT/evGF2.f90
+++ b/src/MBPT/evGF2.f90
@@ -37,19 +37,15 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
   integer                       :: nSCF
   integer                       :: n_diis
   double precision              :: EcBSE(nspin)
-  double precision              :: num
-  double precision              :: eps
   double precision              :: Conv
   double precision              :: rcond
   double precision,allocatable  :: eGF2(:)
   double precision,allocatable  :: eOld(:)
-  double precision,allocatable  :: Sig(:)
+  double precision,allocatable  :: SigC(:)
   double precision,allocatable  :: Z(:)
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: e_diis(:,:)
 
-  integer                       :: i,j,a,b,p
-
 ! Hello world
 
   write(*,*)
@@ -60,7 +56,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
 
 ! Memory allocation
 
-  allocate(Sig(nBas),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
+  allocate(SigC(nBas),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
 
 ! Initialization
 
@@ -80,50 +76,15 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
 
     ! Frequency-dependent second-order contribution
 
-    Sig(:) = 0d0
-    Z(:)   = 0d0
-
-    do p=nC+1,nBas-nR
-      do i=nC+1,nO
-        do j=nC+1,nO
-          do a=nO+1,nBas-nR
-
-            eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
-            num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
-
-            Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
-            Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-          end do
-        end do
-      end do
-    end do
-
-    do p=nC+1,nBas-nR
-      do i=nC+1,nO
-        do a=nO+1,nBas-nR
-          do b=nO+1,nBas-nR
-
-            eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
-            num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
-
-            Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
-            Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-          end do
-        end do
-      end do
-    end do
-
-    Z(:) = 1d0/(1d0 - Z(:))
+    call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
 
     if(linearize) then
 
-      eGF2(:) = eHF(:) + Z(:)*Sig(:)
+      eGF2(:) = eHF(:) + Z(:)*SigC(:)
 
     else
 
-      eGF2(:) = eHF(:) + Sig(:)
+      eGF2(:) = eHF(:) + SigC(:)
 
     end if
 
@@ -131,7 +92,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
 
     ! Print results
 
-    call print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
+    call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF2)
 
     ! DIIS extrapolation
 
diff --git a/src/MBPT/self_energy_GF2.f90 b/src/MBPT/self_energy_GF2.f90
index e4dcdd1..c64d507 100644
--- a/src/MBPT/self_energy_GF2.f90
+++ b/src/MBPT/self_energy_GF2.f90
@@ -1,6 +1,6 @@
-subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
+subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
 
-! Compute GF2 self-energy and its renormalization factor
+! Compute diagonal part of the GF2 self-energy and its renormalization factor
 
   implicit none
   include 'parameters.h'
@@ -16,53 +16,49 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
 ! Local variables
 
   integer                       :: i,j,a,b
-  integer                       :: p,q
+  integer                       :: p
   double precision              :: eps
   double precision              :: num
 
 ! Output variables
 
-  double precision,intent(out)  :: SigC(nBas,nBas)
+  double precision,intent(out)  :: SigC(nBas)
   double precision,intent(out)  :: Z(nBas)
 
 ! Initialize 
 
-  SigC(:,:) = 0d0
-  Z(:)     = 0d0
+  SigC(:) = 0d0
+  Z(:)    = 0d0
 
 ! Compute GF2 self-energy
 
   do p=nC+1,nBas-nR
-    do q=nC+1,nBas-nR
-      do i=nC+1,nO
-        do j=nC+1,nO
-          do a=nO+1,nBas-nR
+    do i=nC+1,nO
+      do j=nC+1,nO
+        do a=nO+1,nBas-nR
 
-            eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
-            num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
+          eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
+          num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
 
-            SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
-            if(p == q) Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+          SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
+          Z(p)    = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
-          end do
         end do
       end do
     end do
   end do
 
   do p=nC+1,nBas-nR
-    do q=nC+1,nBas-nR
-      do i=nC+1,nO
-        do a=nO+1,nBas-nR
-          do b=nO+1,nBas-nR
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
 
-            eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
-            num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
+          eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
+          num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
 
-            SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
-            if(p == q) Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+          SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
+          Z(p)    = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
-          end do
         end do
       end do
     end do
@@ -70,4 +66,4 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
 
   Z(:) = 1d0/(1d0 - Z(:))
 
-end subroutine self_energy_GF2
+end subroutine self_energy_GF2_diag
diff --git a/src/MBPT/self_energy_GF2_diag.f90 b/src/MBPT/self_energy_GF2_diag.f90
new file mode 100644
index 0000000..e4dcdd1
--- /dev/null
+++ b/src/MBPT/self_energy_GF2_diag.f90
@@ -0,0 +1,73 @@
+subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
+
+! Compute GF2 self-energy and its renormalization factor
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eGF2(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+  integer                       :: p,q
+  double precision              :: eps
+  double precision              :: num
+
+! Output variables
+
+  double precision,intent(out)  :: SigC(nBas,nBas)
+  double precision,intent(out)  :: Z(nBas)
+
+! Initialize 
+
+  SigC(:,:) = 0d0
+  Z(:)     = 0d0
+
+! Compute GF2 self-energy
+
+  do p=nC+1,nBas-nR
+    do q=nC+1,nBas-nR
+      do i=nC+1,nO
+        do j=nC+1,nO
+          do a=nO+1,nBas-nR
+
+            eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
+            num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
+
+            SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
+            if(p == q) Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do p=nC+1,nBas-nR
+    do q=nC+1,nBas-nR
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          do b=nO+1,nBas-nR
+
+            eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
+            num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
+
+            SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
+            if(p == q) Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+end subroutine self_energy_GF2