diff --git a/examples/molecule.H2 b/examples/molecule.H2
index 779d849..93d91e8 100644
--- a/examples/molecule.H2
+++ b/examples/molecule.H2
@@ -2,4 +2,4 @@
  2   1     1   0     0
 # Znuc   x            y           z
   H     0.           0.         0.
-  H     0.           0.         1.399
+  H     0.           0.         1.402
diff --git a/examples/molecule.LiH b/examples/molecule.LiH
index 7b1b846..635c08d 100644
--- a/examples/molecule.LiH
+++ b/examples/molecule.LiH
@@ -2,4 +2,4 @@
  2   2     2   0     0
 # Znuc   x            y           z
   Li    0.           0.         0.
-  H     0.           0.         3.5
+  H     0.           0.         3.017
diff --git a/input/basis b/input/basis
index 1ea2746..6f3d2a9 100644
--- a/input/basis
+++ b/input/basis
@@ -1,30 +1,39 @@
-1   6
+1  10
 S   8
-  1   2940.0000000              0.0006800        
-  2    441.2000000              0.0052360        
-  3    100.5000000              0.0266060        
-  4     28.4300000              0.0999930        
-  5      9.1690000              0.2697020        
-  6      3.1960000              0.4514690        
-  7      1.1590000              0.2950740        
-  8      0.1811000              0.0125870        
+  1  24350.0000000              0.0005020        
+  2   3650.0000000              0.0038810        
+  3    829.6000000              0.0199970        
+  4    234.0000000              0.0784180        
+  5     75.6100000              0.2296760        
+  6     26.7300000              0.4327220        
+  7      9.9270000              0.3506420        
+  8      1.1020000             -0.0076450        
 S   8
-  1   2940.0000000             -0.0001230        
-  2    441.2000000             -0.0009660        
-  3    100.5000000             -0.0048310        
-  4     28.4300000             -0.0193140        
-  5      9.1690000             -0.0532800        
-  6      3.1960000             -0.1207230        
-  7      1.1590000             -0.1334350        
-  8      0.1811000              0.5307670        
+  1  24350.0000000             -0.0001180        
+  2   3650.0000000             -0.0009150        
+  3    829.6000000             -0.0047370        
+  4    234.0000000             -0.0192330        
+  5     75.6100000             -0.0603690        
+  6     26.7300000             -0.1425080        
+  7      9.9270000             -0.1777100        
+  8      1.1020000              0.6058360        
 S   1
-  1      0.0589000              1.0000000        
+  1      2.8360000              1.0000000        
+S   1
+  1      0.3782000              1.0000000        
 P   3
-  1      3.6190000              0.0291110        
-  2      0.7110000              0.1693650        
-  3      0.1951000              0.5134580        
+  1     54.7000000              0.0171510        
+  2     12.4300000              0.1076560        
+  3      3.6790000              0.3216810        
 P   1
-  1      0.0601800              1.0000000        
+  1      1.1430000              1.0000000        
+P   1
+  1      0.3300000              1.0000000        
 D   1
-  1      0.2380000              1.0000000
+  1      4.0140000              1.0000000        
+D   1
+  1      1.0960000              1.0000000        
+F   1
+  1      2.5440000              1.0000000
+
 
diff --git a/input/methods b/input/methods
index 2aa67c1..5baa140 100644
--- a/input/methods
+++ b/input/methods
@@ -7,9 +7,9 @@
 # CIS RPA RPAx ppRPA ADC
   F   F   F    F     F
 # GF2 GF3
-  F   F
+  T   F
 # G0W0 evGW qsGW 
-  T    F    F
+  F    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2
diff --git a/input/molecule b/input/molecule
index 6a6f6d1..edeba31 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,4 @@
 # nAt nEla nElb nCore nRyd
-  1   2     2   0     0
+  1   5     5   0     0
 # Znuc   x            y           z
-  Be      0.0          0.0         0.0
+  Ne    0.0          0.0         0.0
diff --git a/input/molecule.xyz b/input/molecule.xyz
index 8023e37..1c70680 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,3 +1,3 @@
   1
 
- Be     0.0000000000    0.0000000000    0.0000000000
+ Ne     0.0000000000    0.0000000000    0.0000000000
diff --git a/input/options b/input/options
index 3be4e0b..5b85bcd 100644
--- a/input/options
+++ b/input/options
@@ -11,6 +11,6 @@
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta     
       256    0.00001    T  5     F      F     T   F   F   F   T    0.000   
 # ACFDT: AC Kx  XBS
-         T  F   T
+         T  F   F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/input/weight b/input/weight
index 1ea2746..6f3d2a9 100644
--- a/input/weight
+++ b/input/weight
@@ -1,30 +1,39 @@
-1   6
+1  10
 S   8
-  1   2940.0000000              0.0006800        
-  2    441.2000000              0.0052360        
-  3    100.5000000              0.0266060        
-  4     28.4300000              0.0999930        
-  5      9.1690000              0.2697020        
-  6      3.1960000              0.4514690        
-  7      1.1590000              0.2950740        
-  8      0.1811000              0.0125870        
+  1  24350.0000000              0.0005020        
+  2   3650.0000000              0.0038810        
+  3    829.6000000              0.0199970        
+  4    234.0000000              0.0784180        
+  5     75.6100000              0.2296760        
+  6     26.7300000              0.4327220        
+  7      9.9270000              0.3506420        
+  8      1.1020000             -0.0076450        
 S   8
-  1   2940.0000000             -0.0001230        
-  2    441.2000000             -0.0009660        
-  3    100.5000000             -0.0048310        
-  4     28.4300000             -0.0193140        
-  5      9.1690000             -0.0532800        
-  6      3.1960000             -0.1207230        
-  7      1.1590000             -0.1334350        
-  8      0.1811000              0.5307670        
+  1  24350.0000000             -0.0001180        
+  2   3650.0000000             -0.0009150        
+  3    829.6000000             -0.0047370        
+  4    234.0000000             -0.0192330        
+  5     75.6100000             -0.0603690        
+  6     26.7300000             -0.1425080        
+  7      9.9270000             -0.1777100        
+  8      1.1020000              0.6058360        
 S   1
-  1      0.0589000              1.0000000        
+  1      2.8360000              1.0000000        
+S   1
+  1      0.3782000              1.0000000        
 P   3
-  1      3.6190000              0.0291110        
-  2      0.7110000              0.1693650        
-  3      0.1951000              0.5134580        
+  1     54.7000000              0.0171510        
+  2     12.4300000              0.1076560        
+  3      3.6790000              0.3216810        
 P   1
-  1      0.0601800              1.0000000        
+  1      1.1430000              1.0000000        
+P   1
+  1      0.3300000              1.0000000        
 D   1
-  1      0.2380000              1.0000000
+  1      4.0140000              1.0000000        
+D   1
+  1      1.0960000              1.0000000        
+F   1
+  1      2.5440000              1.0000000
+
 
diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90
new file mode 100644
index 0000000..f78ef91
--- /dev/null
+++ b/src/QuAcK/G0F2.f90
@@ -0,0 +1,122 @@
+subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
+
+! Perform a one-shot second-order Green function calculation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: linearize
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  double precision,intent(in)   :: e0(nBas)
+  double precision,intent(in)   :: V(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  double precision              :: eps
+  double precision,allocatable  :: eGF2(:)
+  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Z(:)
+
+  integer                       :: i,j,a,b,p
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|     One-shot second-order Green function     |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Memory allocation
+
+  allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas))
+
+
+  ! Frequency-dependent second-order contribution
+
+  Bpp(:,:) = 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) + e0(a) - e0(i) - e0(j)
+
+          Bpp(p,1) = Bpp(p,1) &
+                   + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
+
+        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) + e0(i) - e0(a) - e0(b)
+
+          Bpp(p,2) = Bpp(p,2) &
+                   + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
+
+        end do
+      end do
+    end do
+  end do
+
+  ! Compute the renormalization factor
+
+  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) + e0(a) - e0(i) - e0(j)
+
+          Z(p) = Z(p) - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps**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) + e0(i) - e0(a) - e0(b)
+
+          Z(p) = Z(p) - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps**2
+
+        end do
+      end do
+    end do
+  end do
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+  if(linearize) then
+
+    eGF2(:) = e0(:) + Z(:)*(Bpp(:,1) + Bpp(:,2))
+
+  else
+
+    eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
+
+  end if
+  ! Print results
+
+  call print_G0F2(nBas,nO,e0,eGF2)
+
+end subroutine G0F2
diff --git a/src/QuAcK/G0F3.f90 b/src/QuAcK/G0F3.f90
new file mode 100644
index 0000000..ad972cb
--- /dev/null
+++ b/src/QuAcK/G0F3.f90
@@ -0,0 +1,433 @@
+ subroutine G0F3(renormalization,nBas,nC,nO,nV,nR,V,e0)
+
+! Perform third-order Green function calculation in diagonal approximation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: renormalization
+  integer,intent(in)            :: nBas,nC,nO,nV,nR
+  double precision,intent(in)   :: e0(nBas),V(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  double precision              :: eps,eps1,eps2
+  double precision,allocatable  :: Sig2(:),SigInf(:),Sig3(:),eGF3(:),eOld(:)
+  double precision,allocatable  :: App(:,:),Bpp(:,:),Cpp(:,:),Dpp(:,:)
+  double precision,allocatable  :: Z(:),X2h1p(:),X1h2p(:),Sig2h1p(:),Sig1h2p(:)
+
+  integer                       :: i,j,k,l,a,b,c,d,p
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|    Third-order Green function calculation    |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Memory allocation
+
+  allocate(eGF3(nBas),Sig2(nBas),SigInf(nBas),Sig3(nBas),   &
+           App(nBas,6),Bpp(nBas,2),Cpp(nBas,6),Dpp(nBas,6), &
+           Z(nBas),X2h1p(nBas),X1h2p(nBas),Sig2h1p(nBas),Sig1h2p(nBas))
+
+!------------------------------------------------------------------------
+! Compute third-order frequency-independent contribution
+!------------------------------------------------------------------------
+
+  App(:,:) = 0d0
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+    do k=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+          eps1 = e0(j) + e0(i) - e0(a) - e0(b)
+          eps2 = e0(k) + e0(i) - e0(a) - e0(b)
+
+          App(p,1) = App(p,1) & 
+                   - (2d0*V(p,k,p,j) - V(p,k,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,k,i)/(eps1*eps2)
+
+      enddo
+      enddo
+    enddo
+    enddo
+    enddo
+  enddo
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+      do c=nO+1,nBas-nR
+
+          eps1 = e0(j) + e0(i) - e0(a) - e0(b)
+          eps2 = e0(j) + e0(i) - e0(a) - e0(c)
+
+          App(p,2) = App(p,2) & 
+                   + (2d0*V(p,c,p,b) - V(p,c,b,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(j,i,c,a)/(eps1*eps2)
+
+      enddo
+      enddo
+      enddo
+    enddo
+    enddo
+  enddo
+  
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+      do c=nO+1,nBas-nR
+
+          eps1 = e0(j) + e0(i) - e0(a) - e0(b)
+          eps2 = e0(j)         - e0(c)
+
+          App(p,3) = App(p,3) & 
+                   + (2d0*V(p,c,p,j) - V(p,c,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,c,i)/(eps1*eps2)
+
+      enddo
+      enddo
+      enddo
+    enddo
+    enddo
+  enddo
+
+  App(:,4) = App(:,3)
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+    do k=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+          eps1 = e0(j) + e0(i) - e0(a) - e0(b)
+          eps2 = e0(k)         - e0(b)
+
+          App(p,5) = App(p,5) & 
+                   - (2d0*V(p,b,p,k) - V(p,b,k,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(i,j,k,a)/(eps1*eps2)
+
+      enddo
+      enddo
+    enddo
+    enddo
+    enddo
+  enddo
+  
+  App(:,6) = App(:,5)
+
+! Frequency-independent part of the third-order self-energy
+
+  SigInf(:) = App(:,1) + App(:,2) + App(:,3) + App(:,4) + App(:,5) + App(:,6)
+
+! Frequency-dependent second-order contribution
+
+  Bpp(:,:) = 0d0
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+
+        eps = eGF3(p) + e0(a) - e0(i) - e0(j)
+
+        Bpp(p,1) = Bpp(p,1) & 
+                 + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
+
+      enddo
+    enddo
+    enddo
+  enddo
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+        eps = eGF3(p) + e0(i) - e0(a) - e0(b)
+
+        Bpp(p,2) = Bpp(p,2) &
+                 + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
+
+      enddo
+      enddo
+    enddo
+  enddo
+
+  ! Total second-order Green function
+
+  Sig2(:) = Bpp(:,1) + Bpp(:,2)
+
+  ! Frequency-dependent third-order contribution: "C" terms
+
+  Cpp(:,:) = 0d0
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+      do c=nO+1,nBas-nR
+      do d=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
+          eps2 = eGF3(p) + e0(i) - e0(c) - e0(d)
+
+          Cpp(p,1) = Cpp(p,1) & 
+                   + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,c,d)*V(p,i,c,d)/(eps1*eps2)
+
+      enddo
+      enddo
+      enddo
+      enddo
+    enddo
+  enddo
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+    do k=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
+          eps2 = e0(j)   + e0(k) - e0(a) - e0(b)
+
+          Cpp(p,2) = Cpp(p,2) & 
+                   + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,j,k)*V(p,i,j,k)/(eps1*eps2)
+
+      enddo
+      enddo
+    enddo
+    enddo
+    enddo
+  enddo
+
+  Cpp(:,3) = Cpp(:,2)
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+      do c=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
+          eps2 = e0(i)   + e0(j) - e0(b) - e0(c)
+
+          Cpp(p,4) = Cpp(p,4) & 
+                   + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(i,j,b,c)*V(p,a,b,c)/(eps1*eps2)
+      enddo
+      enddo
+      enddo
+    enddo
+    enddo
+  enddo
+
+  Cpp(:,5) = Cpp(:,4)
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+    do k=nC+1,nO
+    do l=nC+1,nO
+      do a=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
+          eps2 = eGF3(p) + e0(a) - e0(k) - e0(l)
+
+          Cpp(p,6) = Cpp(p,6) & 
+                   - (2d0*V(p,a,k,l) - V(p,a,l,k))*V(k,l,i,j)*V(p,a,i,j)/(eps1*eps2)
+      enddo
+    enddo
+    enddo
+    enddo
+    enddo
+  enddo
+
+  ! Frequency-dependent third-order contribution: "D" terms
+
+  Dpp(:,:) = 0d0
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+      do c=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
+          eps2 = eGF3(p) + e0(j) - e0(b) - e0(c)
+
+          Dpp(p,1) = Dpp(p,1) &
+                   + V(p,i,a,b)*(V(a,j,i,c)*(    V(p,j,c,b) - 2d0*V(p,j,b,c)) &
+                               + V(a,j,c,i)*(    V(p,j,b,c) - 2d0*V(p,j,c,b)))/(eps1*eps2)
+
+          Dpp(p,1) = Dpp(p,1) &
+                   + V(p,i,b,a)*(V(a,j,i,c)*(4d0*V(p,j,b,c) - 2d0*V(p,j,c,b)) &
+                               + V(a,j,c,i)*(    V(p,j,c,b) - 2d0*V(p,j,b,c)))/(eps1*eps2)
+
+      enddo
+      enddo
+      enddo
+    enddo
+    enddo
+  enddo
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+      do c=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(i) - e0(a) - e0(c)
+          eps2 = e0(i)   + e0(j) - e0(a) - e0(b)
+
+          Dpp(p,2) = Dpp(p,2) &
+                   + V(p,i,c,a)*(V(a,b,i,j)*(4d0*V(p,b,c,j) - 2d0*V(p,b,j,c)) &
+                               + V(a,b,j,i)*(    V(p,b,j,c) - 2d0*V(p,b,c,j)))/(eps1*eps2)
+
+          Dpp(p,2) = Dpp(p,2) &
+                   + V(p,i,a,c)*(V(a,b,i,j)*(    V(p,b,j,c) - 2d0*V(p,b,c,j)) &
+                               + V(a,b,j,i)*(    V(p,b,c,j) - 2d0*V(p,b,j,c)))/(eps1*eps2)
+
+      enddo
+      enddo
+      enddo
+    enddo
+    enddo
+  enddo
+
+  Dpp(:,3) = Dpp(:,2)
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+    do k=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(a) - e0(j) - e0(k)
+          eps2 = e0(i)   + e0(j) - e0(a) - e0(b)
+
+          Dpp(p,4) = Dpp(p,4) &
+                   + V(p,a,k,j)*(V(j,i,a,b)*(4d0*V(p,i,k,b) - 2d0*V(p,i,b,k)) &
+                               + V(j,i,b,a)*(    V(p,i,b,k) - 2d0*V(p,i,k,b)))/(eps1*eps2)
+
+          Dpp(p,4) = Dpp(p,4) &
+                   + V(p,a,j,k)*(V(j,i,a,b)*(    V(p,i,b,k) - 2d0*V(p,i,k,b)) &
+                               + V(j,i,b,a)*(    V(p,i,k,b) - 2d0*V(p,i,b,k)))/(eps1*eps2)
+
+      enddo
+      enddo
+    enddo
+    enddo
+    enddo
+  enddo
+
+  Dpp(:,5) = Dpp(:,4)
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+    do j=nC+1,nO
+    do k=nC+1,nO
+      do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+          eps1 = eGF3(p) + e0(a) - e0(i) - e0(k)
+          eps2 = eGF3(p) + e0(b) - e0(j) - e0(k)
+
+          Dpp(p,6) = Dpp(p,6) &
+                   - V(p,a,k,i)*(V(i,b,a,j)*(4d0*V(p,b,k,j) - 2d0*V(p,b,j,k)) &
+                               + V(i,b,j,a)*(    V(p,b,j,k) - 2d0*V(p,b,k,j)))/(eps1*eps2)
+
+          Dpp(p,6) = Dpp(p,6) &
+                   - V(p,a,i,k)*(V(i,b,a,j)*(    V(p,b,j,k) - 2d0*V(p,b,k,j)) &
+                               + V(i,b,j,a)*(    V(p,b,k,j) - 2d0*V(p,b,j,k)))/(eps1*eps2)
+
+      enddo
+      enddo
+    enddo
+    enddo
+    enddo
+  enddo
+
+  ! Compute renormalization factor (if required)
+
+  Z(:) = 1d0
+
+  if(renormalization == 0) then
+
+    Sig3(:) = SigInf(:) &
+            + Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
+            + Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
+
+  elseif(renormalization == 1) then
+
+    Sig3(:) = SigInf(:) &
+            + Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
+            + Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
+
+    Z(:) = Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) & 
+         + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5)
+
+    Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/Sig2(nC+1:nBas-nR)
+    Z(:) = 1d0/(1d0 - Z(:))
+  
+    Sig3(:) = Z(:)*Sig3(:)
+
+  elseif(renormalization == 2) then
+
+    Sig2h1p(:) = Cpp(:,4) + Cpp(:,5) + Cpp(:,6) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
+    Sig1h2p(:) = Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Dpp(:,1) + Dpp(:,2) + Dpp(:,3)
+
+    X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
+    X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
+
+    X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
+    X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
+
+    Sig3(:) = SigInf(:) +                     &
+            + 1d0/(1d0 - X2h1p(:))*Sig2h1p(:) &
+            + 1d0/(1d0 - X1h2p(:))*Sig1h2p(:)
+
+  elseif(renormalization == 3) then
+
+    Sig3(:) = SigInf(:) &
+            + Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
+            + Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
+
+    Sig2h1p(:) = Cpp(:,4) + Cpp(:,5) + Cpp(:,6) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
+    Sig1h2p(:) = Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Dpp(:,1) + Dpp(:,2) + Dpp(:,3)
+
+    X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
+    X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
+
+    X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
+    X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
+
+    Z(:) = X2h1p(:)*Sig2h1p(:) + X1h2p(:)*Sig1h2p(:)
+    Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/(Sig3(nC+1:nBas-nR) - SigInf(nC+1:nBas-nR))
+    Z(:) = 1d0/(1d0 - Z(:))
+
+    Sig3(:) = Z(:)*Sig3(:)
+
+  endif
+
+  ! Total third-order Green function
+
+   eGF3(:) = e0(:) + Sig2(:) + Sig3(:)
+
+  ! Print results
+
+  call print_G0F3(nBas,nO,e0,Z,eGF3)
+
+end subroutine G0F3
diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90
index 7269a5c..7a541d2 100644
--- a/src/QuAcK/G0T0.f90
+++ b/src/QuAcK/G0T0.f90
@@ -77,6 +77,8 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
                           Omega2s(:),X2s(:,:),Y2s(:,:),   & 
                           EcRPA(ispin))
 
+  EcRPA(ispin) = 1d0*EcRPA(ispin)
+
   call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s(:))
   call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s(:))
 
@@ -100,6 +102,8 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
                           Omega2t(:),X2t(:,:),Y2t(:,:),   & 
                           EcRPA(ispin))
 
+  EcRPA(ispin) = 3d0*EcRPA(ispin)
+
   call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:))
   call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:))
 
diff --git a/src/QuAcK/GF2.f90 b/src/QuAcK/GF2.f90
deleted file mode 100644
index 86e8c37..0000000
--- a/src/QuAcK/GF2.f90
+++ /dev/null
@@ -1,137 +0,0 @@
-subroutine GF2(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
-
-! Perform second-order Green function calculation in diagonal approximation
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: maxSCF
-  double precision,intent(in)   :: thresh
-  integer,intent(in)            :: max_diis
-  integer,intent(in)            :: nBas,nC,nO,nV,nR
-  double precision,intent(in)   :: e0(nBas),V(nBas,nBas,nBas,nBas)
-
-! Local variables
-
-  integer                       :: nSCF,n_diis
-  double precision              :: eps
-  double precision              :: Conv
-  double precision              :: rcond
-  double precision,allocatable  :: eGF2(:),eOld(:),Bpp(:,:,:),error_diis(:,:),e_diis(:,:)
-
-  integer                       :: i,j,a,b,p,q
-
-! Hello world
-
-  write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|  Second-order Green function calculation     |'
-  write(*,*)'************************************************'
-  write(*,*)
-
-! Memory allocation
-
-  allocate(Bpp(nBas,nBas,2),eGF2(nBas),eOld(nBas), &
-           error_diis(nBas,max_diis),e_diis(nBas,max_diis))
-
-! Initialization
-
-  Conv            = 1d0
-  nSCF            = 0
-  n_diis          = 0
-  e_diis(:,:)     = 0d0
-  error_diis(:,:) = 0d0
-  eGF2(:)         = e0(:)
-  eOld(:)         = e0(:)
-
-!------------------------------------------------------------------------
-! Main SCF loop
-!------------------------------------------------------------------------
-
-  do while(Conv > thresh .and. nSCF < maxSCF)
-
-    ! Frequency-dependent second-order contribution
-
-    Bpp(:,:,:) = 0d0
-
-    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) + e0(a) - e0(i) - e0(j)
-
-          Bpp(p,q,1) = Bpp(p,q,1) &
-                   + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(q,a,i,j)/eps
-
-        enddo
-      enddo
-      enddo
-    enddo
-    enddo
-
-    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) + e0(i) - e0(a) - e0(b)
-
-          Bpp(p,q,2) = Bpp(p,q,2) &
-                   + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(q,i,a,b)/eps
-
-        enddo
-        enddo
-      enddo
-    enddo
-    enddo
-
-    print*,'Sig2 in GF2'
-    call matout(nBas,nBas,Bpp(:,:,1) + Bpp(:,:,2))
-
-!   eGF2(:) = e0(:) &
-!           + Bpp(:,1) + Bpp(:,2)
-
-    Conv = maxval(abs(eGF2 - eOld))
-
-    ! DIIS extrapolation
-
-    n_diis = min(n_diis+1,max_diis)
-    call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2)
-
-!    Reset DIIS if required
-
-    if(abs(rcond) < 1d-15) n_diis = 0
-
-    eOld = eGF2
-
-    ! Print results
-
-    call print_GF2(nBas,nO,nSCF,Conv,e0,eGF2)
-
-    ! Increment
-
-    nSCF = nSCF + 1
-
-  enddo
-!------------------------------------------------------------------------
-! End main SCF loop
-!------------------------------------------------------------------------
-
-! Did it actually converge?
-
-  if(nSCF == maxSCF+1) then
-
-    write(*,*)
-    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
-    write(*,*)'                 Convergence failed                 '
-    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
-    write(*,*)
-
-  endif
-
-end subroutine GF2
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 05e6384..5d6d9ea 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -10,7 +10,7 @@ program QuAcK
   logical                       :: do_ring_CCD,do_ladder_CCD
   logical                       :: doCIS,doRPA,doRPAx
   logical                       :: doppRPA,doADC
-  logical                       :: doGF2,doGF3
+  logical                       :: doG0F2,doevGF2,doG0F3,doevGF3
   logical                       :: doG0W0,doevGW,doqsGW
   logical                       :: doG0T0,doevGT,doqsGT
   logical                       :: doMCMP2,doMinMCMP2
@@ -114,15 +114,15 @@ program QuAcK
 
 ! Which calculations do you want to do?
 
-  call read_methods(doRHF,doUHF,doMOM,         &
-                    doMP2,doMP3,doMP2F12,      &
-                    doCCD,doCCSD,doCCSDT,      &
-                    do_ring_CCD,do_ladder_CCD, &
-                    doCIS,doRPA,doRPAx,        & 
-                    doppRPA,doADC,             &
-                    doGF2,doGF3,               &
-                    doG0W0,doevGW,doqsGW,      &
-                    doG0T0,doevGT,doqsGT,      &
+  call read_methods(doRHF,doUHF,doMOM,             &
+                    doMP2,doMP3,doMP2F12,          &
+                    doCCD,doCCSD,doCCSDT,          &
+                    do_ring_CCD,do_ladder_CCD,     &
+                    doCIS,doRPA,doRPAx,            & 
+                    doppRPA,doADC,                 &
+                    doG0F2,doevGF2,doG0F3,doevGF3, &
+                    doG0W0,doevGW,doqsGW,          &
+                    doG0T0,doevGT,doqsGT,          &
                     doMCMP2)
 
 ! Read options for methods
@@ -507,13 +507,13 @@ program QuAcK
   end if
 
 !------------------------------------------------------------------------
-! Compute GF2 electronic binding energies
+! Compute G0F2 electronic binding energies
 !------------------------------------------------------------------------
 
-  if(doGF2) then
+  if(doG0F2) then
 
     call cpu_time(start_GF2)
-    call GF2_diag(maxSCF_GF,thresh_GF,n_diis_GF,linearize,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
+    call G0F2(maxSCF_GF,thresh_GF,n_diis_GF,linearize,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
     call cpu_time(end_GF2)
 
     t_GF2 = end_GF2 - start_GF2
@@ -523,13 +523,45 @@ program QuAcK
   end if
 
 !------------------------------------------------------------------------
-! Compute GF3 electronic binding energies
+! Compute evGF2 electronic binding energies
 !------------------------------------------------------------------------
 
-  if(doGF3) then
+  if(doevGF2) then
+
+    call cpu_time(start_GF2)
+    call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linearize,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
+    call cpu_time(end_GF2)
+
+    t_GF2 = end_GF2 - start_GF2
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Compute G0F3 electronic binding energies
+!------------------------------------------------------------------------
+
+  if(doG0F3) then
 
     call cpu_time(start_GF3)
-    call GF3_diag(maxSCF_GF,thresh_GF,n_diis_GF,renormalization,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
+    call G0F3(renormalization,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
+    call cpu_time(end_GF3)
+
+    t_GF3 = end_GF3 - start_GF3
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Compute evGF3 electronic binding energies
+!------------------------------------------------------------------------
+
+  if(doevGF3) then
+
+    call cpu_time(start_GF3)
+    call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormalization,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
     call cpu_time(end_GF3)
 
     t_GF3 = end_GF3 - start_GF3
diff --git a/src/QuAcK/GF2_diag.f90 b/src/QuAcK/evGF2.f90
similarity index 94%
rename from src/QuAcK/GF2_diag.f90
rename to src/QuAcK/evGF2.f90
index 28ddb63..e0abc77 100644
--- a/src/QuAcK/GF2_diag.f90
+++ b/src/QuAcK/evGF2.f90
@@ -1,6 +1,6 @@
-subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
+subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
-! Perform second-order Green function calculation in diagonal approximation
+! Perform eigenvalue self-consistent second-order Green function calculation
 
   implicit none
   include 'parameters.h'
@@ -145,7 +145,7 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
     ! Print results
 
-    call print_GF2(nBas,nO,nSCF,Conv,e0,eGF2)
+    call print_evGF2(nBas,nO,nSCF,Conv,e0,eGF2)
 
     ! DIIS extrapolation
 
@@ -177,4 +177,4 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
   end if
 
-end subroutine GF2_diag
+end subroutine evGF2
diff --git a/src/QuAcK/GF3_diag.f90 b/src/QuAcK/evGF3.f90
similarity index 98%
rename from src/QuAcK/GF3_diag.f90
rename to src/QuAcK/evGF3.f90
index 67ad523..b9b6588 100644
--- a/src/QuAcK/GF3_diag.f90
+++ b/src/QuAcK/evGF3.f90
@@ -1,4 +1,4 @@
- subroutine GF3_diag(maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
+ subroutine evGF3(maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
 
 ! Perform third-order Green function calculation in diagonal approximation
 
@@ -455,7 +455,7 @@
 
     ! Print results
 
-    call print_GF3(nBas,nO,nSCF,Conv,e0,Z,eGF3)
+    call print_evGF3(nBas,nO,nSCF,Conv,e0,Z,eGF3)
 
     ! DIIS extrapolation
 
@@ -489,4 +489,4 @@
 
   endif
 
-end subroutine GF3_diag
+end subroutine evGF3
diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90
index df0f510..b9dfdb6 100644
--- a/src/QuAcK/excitation_density_Tmatrix.f90
+++ b/src/QuAcK/excitation_density_Tmatrix.f90
@@ -43,7 +43,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
           cd = 0
           do c=nO+1,nBas-nR
-           do d=nO+1,c
+           do d=c,nBas-nR
               cd = cd + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
 !                          + ERI(p,i,c,d)*X1(cd,ab) 
@@ -54,7 +54,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
           kl = 0
           do k=nC+1,nO
-           do l=nC+1,k
+           do l=k,nO
               kl = kl + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
 !                          + ERI(p,i,k,l)*Y1(kl,ab)
@@ -71,7 +71,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
           cd = 0
           do c=nO+1,nBas-nR
-           do d=nO+1,c
+           do d=c,nBas-nR
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
 !                          + ERI(p,nO+a,c,d)*X2(cd,ij)
@@ -83,7 +83,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
           kl = 0
           do k=nC+1,nO
-           do l=nC+1,k
+           do l=k,nO
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
 !                          + ERI(p,nO+a,k,l)*Y2(kl,ij)  
@@ -105,65 +105,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
   if(ispin == 2) then
 
-    do p=nC+1,nBas-nR
-
-      do i=nC+1,nO
-        do ab=1,nVV
-
-          cd = 0
-          do c=nO+1,nBas-nR
-           do d=nO+1,c-1
-              cd = cd + 1
-              rho1(p,i,ab) = rho1(p,i,ab) & 
-                           + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
-            end do
-          end do
-
-          kl = 0
-          do k=nC+1,nO
-           do l=nC+1,k-1
-              kl = kl + 1
-              rho1(p,i,ab) = rho1(p,i,ab) & 
-                           + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
-            end do
-          end do
-
-        end do
-      end do
-
-      do a=1,nV-nR
-        do ij=1,nOO
-
-          cd = 0
-          do c=nO+1,nBas-nR
-           do d=nO+1,c-1
-              cd = cd + 1
-              rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
-            end do
-          end do
-
-          kl = 0
-          do k=nC+1,nO
-           do l=nC+1,k-1
-              kl = kl + 1
-              rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
-            end do
-          end do
-
-        end do
-      end do
-    end do
-
-  end if
-
-!----------------------------------------------
-! Spinorbital basis
-!----------------------------------------------
-
-  if(ispin == 3) then
-
     do p=nC+1,nBas-nR
 
       do i=nC+1,nO
@@ -195,7 +136,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
           cd = 0
           do c=nO+1,nBas-nR
-            do d=c+1,nBas-nR
+           do d=c+1,nBas-nR
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
                            + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
@@ -213,6 +154,65 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
         end do
       end do
+    end do
+
+  end if
+
+!----------------------------------------------
+! Spinorbital basis
+!----------------------------------------------
+
+  if(ispin == 3) then
+
+    do p=nC+1,nBas-nR
+
+      do i=nC+1,nO
+        do ab=1,nVV
+
+          cd = 0
+          do c=nO+1,nBas-nR
+           do d=c+1,nBas-nR
+              cd = cd + 1
+              rho1(p,i,ab)  = rho1(p,i,ab) & 
+                            + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
+            end do
+          end do
+
+          kl = 0
+          do k=nC+1,nO
+           do l=k+1,nO
+              kl = kl + 1
+              rho1(p,i,ab)  = rho1(p,i,ab) & 
+                            + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
+            end do
+          end do
+
+        end do
+      end do
+
+      do a=1,nV-nR
+        do ij=1,nOO
+
+          cd = 0
+          do c=nO+1,nBas-nR
+            do d=c+1,nBas-nR
+              cd = cd + 1
+              rho2(p,a,ij)  = rho2(p,a,ij) & 
+                            + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
+            end do
+          end do
+
+          kl = 0
+          do k=nC+1,nO
+           do l=k+1,nO
+              kl = kl + 1
+              rho2(p,a,ij)  = rho2(p,a,ij) & 
+                            + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
+            end do
+          end do
+
+        end do
+      end do
 
     end do
 
diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90
index e9d18a5..baa2079 100644
--- a/src/QuAcK/linear_response_B_pp.f90
+++ b/src/QuAcK/linear_response_B_pp.f90
@@ -26,11 +26,11 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 
     ab = 0
     do a=nO+1,nBas-nR
-     do b=nO+1,a
+     do b=a,nBas-nR
         ab = ab + 1
         ij = 0
         do i=nC+1,nO
-         do j=nC+1,i
+         do j=i,nO
             ij = ij + 1
  
             B_pp(ab,ij) = (ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
@@ -48,11 +48,11 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 
     ab = 0
     do a=nO+1,nBas-nR
-     do b=nO+1,a-1
+     do b=a+1,nBas-nR
         ab = ab + 1
         ij = 0
         do i=nC+1,nO
-         do j=nC+1,i-1
+         do j=i+1,nO
             ij = ij + 1
  
             B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i)
diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90
index 72e56f2..b82f1a3 100644
--- a/src/QuAcK/linear_response_C_pp.f90
+++ b/src/QuAcK/linear_response_C_pp.f90
@@ -32,11 +32,11 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 
     ab = 0
     do a=nO+1,nBas-nR
-     do b=nO+1,a
+     do b=a,nBas-nR
         ab = ab + 1
         cd = 0
         do c=nO+1,nBas-nR
-         do d=nO+1,c
+         do d=c,nBas-nR
             cd = cd + 1
  
             C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
@@ -55,11 +55,11 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 
     ab = 0
     do a=nO+1,nBas-nR
-     do b=nO+1,a-1
+     do b=a+1,nBas-nR
         ab = ab + 1
         cd = 0
         do c=nO+1,nBas-nR
-         do d=nO+1,c-1
+         do d=c+1,nBas-nR
             cd = cd + 1
  
             C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90
index fec3d10..1b38222 100644
--- a/src/QuAcK/linear_response_D_pp.f90
+++ b/src/QuAcK/linear_response_D_pp.f90
@@ -32,11 +32,11 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 
     ij = 0
     do i=nC+1,nO
-     do j=nC+1,i
+     do j=i,nO
         ij = ij + 1
         kl = 0
         do k=nC+1,nO
-         do l=nC+1,k
+         do l=k,nO
             kl = kl + 1
  
             D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
@@ -55,11 +55,11 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 
     ij = 0
     do i=nC+1,nO
-     do j=nC+1,i-1
+     do j=i+1,nO
         ij = ij + 1
         kl = 0
         do k=nC+1,nO
-         do l=nC+1,k-1
+         do l=k+1,nO
             kl = kl + 1
  
             D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
diff --git a/src/QuAcK/print_G0F2.f90 b/src/QuAcK/print_G0F2.f90
new file mode 100644
index 0000000..479b959
--- /dev/null
+++ b/src/QuAcK/print_G0F2.f90
@@ -0,0 +1,45 @@
+subroutine print_G0F2(nBas,nO,eHF,eGF2)
+
+! Print one-electron energies and other stuff for G0F2
+
+  implicit none
+  include 'parameters.h'
+
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: eGF2(nBas)
+
+  integer                            :: x
+  integer                            :: HOMO
+  integer                            :: LUMO
+  double precision                   :: Gap
+
+! HOMO and LUMO
+
+  HOMO = nO
+  LUMO = HOMO + 1
+  Gap = eGF2(LUMO)-eGF2(HOMO)
+
+! Dump results
+
+  write(*,*)'-------------------------------------------'
+  write(*,*)' Frequency-dependent G0F2 calculation'
+  write(*,*)'-------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','e_G0F2 (eV)','|'
+  write(*,*)'-------------------------------------------'
+
+  do x=1,nBas
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+    '|',x,'|',eHF(x)*HaToeV,'|',eGF2(x)*HaToeV,'|'
+  enddo
+
+  write(*,*)'-------------------------------------------'
+  write(*,'(2X,A27,F15.6)') 'G0F2  HOMO      energy (eV):',eGF2(HOMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'G0F2  LUMO      energy (eV):',eGF2(LUMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'G0F2  HOMO-LUMO gap    (eV):',Gap*HaToeV
+  write(*,*)'-------------------------------------------'
+  write(*,*)
+
+end subroutine print_G0F2
diff --git a/src/QuAcK/print_G0F3.f90 b/src/QuAcK/print_G0F3.f90
new file mode 100644
index 0000000..09da736
--- /dev/null
+++ b/src/QuAcK/print_G0F3.f90
@@ -0,0 +1,41 @@
+subroutine print_G0F3(nBas,nO,eHF,Z,eGF3)
+
+! Print one-electron energies and other stuff for GF3
+
+  implicit none
+  include 'parameters.h'
+
+  integer,intent(in)                 :: nBas,nO
+  double precision,intent(in)        :: eHF(nBas),eGF3(nBas),Z(nBas)
+
+  integer                            :: x,HOMO,LUMO
+  double precision                   :: Gap
+
+! HOMO and LUMO
+
+  HOMO = nO
+  LUMO = HOMO + 1
+  Gap = eGF3(LUMO)-eGF3(HOMO)
+
+! Dump results
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)' Frequency-dependent G0F3 calculation'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','Z','|','e_G0F3 (eV)','|'
+  write(*,*)'-------------------------------------------------------------'
+
+  do x=1,nBas
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+    '|',x,'|',eHF(x)*HaToeV,'|',Z(x),'|',eGF3(x)*HaToeV,'|'
+  enddo
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(2X,A27,F15.6)') 'G0F3  HOMO      energy (eV):',eGF3(HOMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'G0F3  LUMO      energy (eV):',eGF3(LUMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'G0F3  HOMO-LUMO gap    (eV):',Gap*HaToeV
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+end subroutine print_G0F3
diff --git a/src/QuAcK/print_GF2.f90 b/src/QuAcK/print_evGF2.f90
similarity index 72%
rename from src/QuAcK/print_GF2.f90
rename to src/QuAcK/print_evGF2.f90
index 98a3933..050dc7c 100644
--- a/src/QuAcK/print_GF2.f90
+++ b/src/QuAcK/print_evGF2.f90
@@ -1,4 +1,4 @@
-subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
+subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2)
 
 ! Print one-electron energies and other stuff for GF2
 
@@ -20,10 +20,10 @@ subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
 ! Dump results
 
   write(*,*)'-------------------------------------------'
-  write(*,*)' Frequency-dependent diagonal GF2 calculation'
+  write(*,*)' Frequency-dependent evGF2 calculation'
   write(*,*)'-------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
-            '|','#','|','e_HF (eV)','|','e_GF2 (eV)','|'
+            '|','#','|','e_HF (eV)','|','e_evGF2 (eV)','|'
   write(*,*)'-------------------------------------------'
 
   do x=1,nBas
@@ -35,10 +35,10 @@ subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
   write(*,'(2X,A10,I3)')   'Iteration ',nSCF
   write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
   write(*,*)'-------------------------------------------'
-  write(*,'(2X,A27,F15.6)') 'GF2  HOMO      energy (eV):',eGF2(HOMO)*HaToeV
-  write(*,'(2X,A27,F15.6)') 'GF2  LUMO      energy (eV):',eGF2(LUMO)*HaToeV
-  write(*,'(2X,A27,F15.6)') 'GF2  HOMO-LUMO gap    (eV):',Gap*HaToeV
+  write(*,'(2X,A27,F15.6)') 'evGF2  HOMO      energy (eV):',eGF2(HOMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'evGF2  LUMO      energy (eV):',eGF2(LUMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'evGF2  HOMO-LUMO gap    (eV):',Gap*HaToeV
   write(*,*)'-------------------------------------------'
   write(*,*)
 
-end subroutine print_GF2
+end subroutine print_evGF2
diff --git a/src/QuAcK/print_GF3.f90 b/src/QuAcK/print_evGF3.f90
similarity index 73%
rename from src/QuAcK/print_GF3.f90
rename to src/QuAcK/print_evGF3.f90
index cf0634b..3e7f740 100644
--- a/src/QuAcK/print_GF3.f90
+++ b/src/QuAcK/print_evGF3.f90
@@ -1,4 +1,4 @@
-subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
+subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
 
 ! Print one-electron energies and other stuff for GF3
 
@@ -20,10 +20,10 @@ subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
 ! Dump results
 
   write(*,*)'-------------------------------------------------------------'
-  write(*,*)' Frequency-dependent diagonal GF3 calculation'
+  write(*,*)' Frequency-dependent diagonal evGF3 calculation'
   write(*,*)'-------------------------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
-            '|','#','|','e_HF (eV)','|','Z','|','e_GF3 (eV)','|'
+            '|','#','|','e_HF (eV)','|','Z','|','e_evGF3 (eV)','|'
   write(*,*)'-------------------------------------------------------------'
 
   do x=1,nBas
@@ -35,10 +35,10 @@ subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
   write(*,'(2X,A10,I3)')   'Iteration ',nSCF
   write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
   write(*,*)'-------------------------------------------------------------'
-  write(*,'(2X,A27,F15.6)') 'GF3  HOMO      energy (eV):',eGF3(HOMO)*HaToeV
-  write(*,'(2X,A27,F15.6)') 'GF3  LUMO      energy (eV):',eGF3(LUMO)*HaToeV
-  write(*,'(2X,A27,F15.6)') 'GF3  HOMO-LUMO gap    (eV):',Gap*HaToeV
+  write(*,'(2X,A27,F15.6)') 'evGF3  HOMO      energy (eV):',eGF3(HOMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'evGF3  LUMO      energy (eV):',eGF3(LUMO)*HaToeV
+  write(*,'(2X,A27,F15.6)') 'evGF3  HOMO-LUMO gap    (eV):',Gap*HaToeV
   write(*,*)'-------------------------------------------------------------'
   write(*,*)
 
-end subroutine print_GF3
+end subroutine print_evGF3
diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90
index df6991f..e8e89fa 100644
--- a/src/QuAcK/read_methods.f90
+++ b/src/QuAcK/read_methods.f90
@@ -1,12 +1,12 @@
-subroutine read_methods(doRHF,doUHF,doMOM,         & 
-                        doMP2,doMP3,doMP2F12,      & 
-                        doCCD,doCCSD,doCCSDT,      & 
-                        do_ring_CCD,do_ladder_CCD, &
-                        doCIS,doRPA,doRPAx,        & 
-                        doppRPA,doADC,             & 
-                        doGF2,doGF3,               & 
-                        doG0W0,doevGW,doqsGW,      & 
-                        doG0T0,doevGT,doqsGT,      & 
+subroutine read_methods(doRHF,doUHF,doMOM,             & 
+                        doMP2,doMP3,doMP2F12,          & 
+                        doCCD,doCCSD,doCCSDT,          & 
+                        do_ring_CCD,do_ladder_CCD,     &
+                        doCIS,doRPA,doRPAx,            & 
+                        doppRPA,doADC,                 & 
+                        doG0F2,doevGF2,doG0F3,doevGF3, & 
+                        doG0W0,doevGW,doqsGW,          & 
+                        doG0T0,doevGT,doqsGT,          & 
                         doMCMP2)
 
 ! Read desired methods 
@@ -20,7 +20,7 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
   logical,intent(out)           :: doCCD,doCCSD,doCCSDT
   logical,intent(out)           :: do_ring_CCD,do_ladder_CCD
   logical,intent(out)           :: doCIS,doRPA,doRPAx,doppRPA,doADC
-  logical,intent(out)           :: doGF2,doGF3  
+  logical,intent(out)           :: doG0F2,doevGF2,doG0F3,doevGF3  
   logical,intent(out)           :: doG0W0,doevGW,doqsGW
   logical,intent(out)           :: doG0T0,doevGT,doqsGT
   logical,intent(out)           :: doMCMP2
@@ -56,8 +56,10 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
   doppRPA = .false.
   doADC   = .false.
 
-  doGF2 = .false.
-  doGF3 = .false.
+  doG0F2  = .false.
+  doevGF2 = .false.
+  doG0F3  = .false.
+  doevGF3 = .false.
 
   doG0W0 = .false.
   doevGT = .false.
@@ -108,9 +110,11 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
 ! Read Green function methods
 
   read(1,*) 
-  read(1,*) answer1,answer2
-  if(answer1 == 'T') doGF2 = .true.
-  if(answer2 == 'T') doGF3 = .true.
+  read(1,*) answer1,answer2,answer3,answer4
+  if(answer1 == 'T') doG0F2  = .true.
+  if(answer2 == 'T') doevGF2 = .true.
+  if(answer3 == 'T') doG0F3  = .true.
+  if(answer4 == 'T') doevGF3 = .true.
 
 ! Read GW methods
 
diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90
index 14ff12a..1eb76fe 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix.f90
@@ -42,7 +42,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do i=nC+1,nO
       do cd=1,nVVs
         eps = e(p) + e(i) - Omega1s(cd)
-        Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2
+        Z(p) = Z(p) + 1d0*(rho1s(p,i,cd)/eps)**2
       enddo
     enddo
   enddo
@@ -53,7 +53,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do a=1,nV-nR
       do kl=1,nOOs
         eps = e(p) + e(nO+a) - Omega2s(kl)
-        Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2
+        Z(p) = Z(p) + 1d0*(rho2s(p,a,kl)/eps)**2
       enddo
     enddo
   enddo
@@ -68,7 +68,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do i=nC+1,nO
       do cd=1,nVVt
         eps = e(p) + e(i) - Omega1t(cd)
-        Z(p) = Z(p) - 2d0*rho1t(p,i,cd)**2/eps**2
+        Z(p) = Z(p) + 1d0*(rho1t(p,i,cd)/eps)**2
       enddo
     enddo
   enddo
@@ -79,13 +79,13 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do a=1,nV-nR
       do kl=1,nOOt
         eps = e(p) + e(nO+a) - Omega2t(kl)
-        Z(p) = Z(p) - 2d0*rho2t(p,a,kl)**2/eps**2
+        Z(p) = Z(p) + 1d0*(rho2t(p,a,kl)/eps)**2
       enddo
     enddo
   enddo
 
 ! Compute renormalization factor from derivative of SigT
  
-  Z(:) = 1d0/(1d0 - Z(:))
+  Z(:) = 1d0/(1d0 + Z(:))
 
 end subroutine renormalization_factor_Tmatrix
diff --git a/src/QuAcK/renormalization_factor_Tmatrix_so.f90 b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
index aa3692a..c545ba2 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix_so.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
@@ -39,12 +39,8 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
   do p=nC+1,nBas-nR
     do i=nC+1,nO
       do cd=1,nVV
-!     do c=nO+1,nBas-nR
-!       do d=c+1,nBas-nR
-!         cd = cd + 1
-          eps = e(p) + e(i) - Omega1(cd)
-          Z(p) = Z(p) - rho1(p,i,cd)**2/eps**2
-!       enddo
+        eps  = e(p) + e(i) - Omega1(cd)
+        Z(p) = Z(p) + (rho1(p,i,cd)/eps)**2
       enddo
     enddo
   enddo
@@ -54,18 +50,14 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
   do p=nC+1,nBas-nR
     do a=1,nV-nR
       do kl=1,nOO
-!     do k=nC+1,nO
-!       do l=k+1,nO
-!         kl = kl + 1
-          eps = e(p) + e(nO+a) - Omega2(kl)
-          Z(p) = Z(p) - rho2(p,a,kl)**2/eps**2
-!       enddo
+        eps  = e(p) + e(nO+a) - Omega2(kl)
+        Z(p) = Z(p) + (rho2(p,a,kl)/eps)**2
       enddo
     enddo
   enddo
 
 ! Compute renormalization factor from derivative of SigT
  
-  Z(:) = 1d0/(1d0 - Z(:))
+  Z(:) = 1d0/(1d0 + Z(:))
 
 end subroutine renormalization_factor_Tmatrix_so
diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90
index b793985..72d7963 100644
--- a/src/QuAcK/self_energy_Tmatrix_diag.f90
+++ b/src/QuAcK/self_energy_Tmatrix_diag.f90
@@ -46,7 +46,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do i=nC+1,nO
       do cd=1,nVVs
         eps = e(p) + e(i) - Omega1s(cd)
-        SigT(p) = SigT(p) + 2d0*rho1s(p,i,cd)**2/eps
+        SigT(p) = SigT(p) + 1d0*rho1s(p,i,cd)**2/eps
       enddo
     enddo
   enddo
@@ -57,7 +57,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do a=1,nV-nR
       do kl=1,nOOs
         eps = e(p) + e(nO+a) - Omega2s(kl)
-        SigT(p) = SigT(p) + 2d0*rho2s(p,a,kl)**2/eps
+        SigT(p) = SigT(p) + 1d0*rho2s(p,a,kl)**2/eps
       enddo
     enddo
   enddo
@@ -72,7 +72,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do i=nC+1,nO
       do cd=1,nVVt
         eps = e(p) + e(i) - Omega1t(cd)
-        SigT(p) = SigT(p) + 2d0*rho1t(p,i,cd)**2/eps
+        SigT(p) = SigT(p) + 1d0*rho1t(p,i,cd)**2/eps
       enddo
     enddo
   enddo
@@ -83,7 +83,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do a=1,nV-nR
       do kl=1,nOOt
         eps = e(p) + e(nO+a) - Omega2t(kl)
-        SigT(p) = SigT(p) + 2d0*rho2t(p,a,kl)**2/eps
+        SigT(p) = SigT(p) + 1d0*rho2t(p,a,kl)**2/eps
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/self_energy_Tmatrix_diag_so.f90 b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
index ab96502..edef028 100644
--- a/src/QuAcK/self_energy_Tmatrix_diag_so.f90
+++ b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
@@ -43,12 +43,8 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
   do p=nC+1,nBas-nR
     do i=nC+1,nO
       do cd=1,nVV
-!     do c=nO+1,nBas-nR
-!       do d=c+1,nBas-nR
-!         cd = cd + 1
-          eps = e(p) + e(i) - Omega1(cd)
-          SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
-!       enddo
+        eps     = e(p)    + e(i) - Omega1(cd)
+        SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
       enddo
     enddo
   enddo
@@ -58,12 +54,8 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
   do p=nC+1,nBas-nR
     do a=1,nV-nR
       do kl=1,nOO
-!     do k=nC+1,nO
-!       do l=k+1,nO
-!         kl = kl + 1
-          eps = e(p) + e(nO+a) - Omega2(kl)
-          SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
-!       enddo
+        eps     = e(p)    + e(nO+a) - Omega2(kl)
+        SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90
index 4e829af..03b1ba9 100644
--- a/src/QuAcK/soG0T0.f90
+++ b/src/QuAcK/soG0T0.f90
@@ -89,7 +89,7 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 ! Compute excitation densities for the T-matrix
 
   call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & 
-                                  X1(:,:),Y1(:,:),rho1(:,:,:),             & 
+                                  X1(:,:),Y1(:,:),rho1(:,:,:), &
                                   X2(:,:),Y2(:,:),rho2(:,:,:))
 
 !----------------------------------------------
@@ -97,13 +97,15 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 !----------------------------------------------
 
   call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
-                                Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & 
+                                Omega1(:),rho1(:,:,:), &
+                                Omega2(:),rho2(:,:,:), &
                                 SigT(:))
 
 ! Compute renormalization factor for T-matrix self-energy
 
   call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
-                                      Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & 
+                                      Omega1(:),rho1(:,:,:), &
+                                      Omega2(:),rho2(:,:,:), &
                                       Z(:))
 
 !----------------------------------------------