diff --git a/examples/molecule.F2 b/examples/molecule.F2
index 6f3bd1d..7ef2def 100644
--- a/examples/molecule.F2
+++ b/examples/molecule.F2
@@ -2,4 +2,4 @@
   2   9   9  0      0
 # Znuc   x            y           z
   F     0.           0.          0.000
-  F     0.           0.          2.668
+  F     0.           0.          2.68454493
diff --git a/examples/molecule.LiH b/examples/molecule.LiH
index 771beea..c38d7a8 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.061356
+  H     0.           0.         3.09839495
diff --git a/input/basis b/input/basis
index 27724b4..11cb61a 100644
--- a/input/basis
+++ b/input/basis
@@ -1,37 +1,76 @@
-1  10 
-S   8   1.00
-  24350.0000000              0.0005020        
-   3650.0000000              0.0038810        
-    829.6000000              0.0199970        
-    234.0000000              0.0784180        
-     75.6100000              0.2296760        
-     26.7300000              0.4327220        
-      9.9270000              0.3506420        
-      1.1020000             -0.0076450        
-S   8   1.00
-  24350.0000000             -0.0001180        
-   3650.0000000             -0.0009150        
-    829.6000000             -0.0047370        
-    234.0000000             -0.0192330        
-     75.6100000             -0.0603690        
-     26.7300000             -0.1425080        
-      9.9270000             -0.1777100        
-      1.1020000              0.6058360        
+1     9
+S   9   1.00
+      1.471000D+04           7.210000D-04
+      2.207000D+03           5.553000D-03
+      5.028000D+02           2.826700D-02
+      1.426000D+02           1.064440D-01
+      4.647000D+01           2.868140D-01
+      1.670000D+01           4.486410D-01
+      6.356000D+00           2.647610D-01
+      1.316000D+00           1.533300D-02
+      3.897000D-01          -2.332000D-03
+S   9   1.00
+      1.471000D+04          -1.650000D-04
+      2.207000D+03          -1.308000D-03
+      5.028000D+02          -6.495000D-03
+      1.426000D+02          -2.669100D-02
+      4.647000D+01          -7.369000D-02
+      1.670000D+01          -1.707760D-01
+      6.356000D+00          -1.123270D-01
+      1.316000D+00           5.628140D-01
+      3.897000D-01           5.687780D-01
 S   1   1.00
-      2.8360000              1.0000000        
+      3.897000D-01           1.000000D+00
 S   1   1.00
-      0.3782000              1.0000000        
-P   3   1.00
-     54.7000000              0.0171510        
-     12.4300000              0.1076560        
-      3.6790000              0.3216810        
+      0.0986300              1.0000000
+P   4   1.00
+      2.267000D+01           4.487800D-02
+      4.977000D+00           2.357180D-01
+      1.347000D+00           5.085210D-01
+      3.471000D-01           4.581200D-01
 P   1   1.00
-      1.1430000              1.0000000        
+      3.471000D-01           1.000000D+00
 P   1   1.00
-      0.3300000              1.0000000        
+      0.0850200              1.0000000
 D   1   1.00
-      4.0140000              1.0000000        
+      1.640000D+00           1.0000000
 D   1   1.00
-      1.0960000              1.0000000        
-F   1   1.00
-      2.5440000              1.0000000   
+      0.4640000              1.0000000
+2   9
+S   9   1.00
+      1.471000D+04           7.210000D-04
+      2.207000D+03           5.553000D-03
+      5.028000D+02           2.826700D-02
+      1.426000D+02           1.064440D-01
+      4.647000D+01           2.868140D-01
+      1.670000D+01           4.486410D-01
+      6.356000D+00           2.647610D-01
+      1.316000D+00           1.533300D-02
+      3.897000D-01          -2.332000D-03
+S   9   1.00
+      1.471000D+04          -1.650000D-04
+      2.207000D+03          -1.308000D-03
+      5.028000D+02          -6.495000D-03
+      1.426000D+02          -2.669100D-02
+      4.647000D+01          -7.369000D-02
+      1.670000D+01          -1.707760D-01
+      6.356000D+00          -1.123270D-01
+      1.316000D+00           5.628140D-01
+      3.897000D-01           5.687780D-01
+S   1   1.00
+      3.897000D-01           1.000000D+00
+S   1   1.00
+      0.0986300              1.0000000
+P   4   1.00
+      2.267000D+01           4.487800D-02
+      4.977000D+00           2.357180D-01
+      1.347000D+00           5.085210D-01
+      3.471000D-01           4.581200D-01
+P   1   1.00
+      3.471000D-01           1.000000D+00
+P   1   1.00
+      0.0850200              1.0000000
+D   1   1.00
+      1.640000D+00           1.0000000
+D   1   1.00
+      0.4640000              1.0000000
diff --git a/input/molecule b/input/molecule
index edeba31..7ef2def 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,5 @@
-# nAt nEla nElb nCore nRyd
-  1   5     5   0     0
+# nAt nEl nCore nRyd
+  2   9   9  0      0
 # Znuc   x            y           z
-  Ne    0.0          0.0         0.0
+  F     0.           0.          0.000
+  F     0.           0.          2.68454493
diff --git a/input/options b/input/options
index bc5c578..5d37d6a 100644
--- a/input/options
+++ b/input/options
@@ -9,6 +9,6 @@
 # GF: maxSCF thresh  DIIS n_diis renormalization
       64     0.00001 T    10      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta
-      256    0.0000001  T  5      F      F     T   F   F   F   F         0.000
+      256    0.0000001  T  5      F      F     T   F   F   F   T         0.000
 # 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 27724b4..11cb61a 100644
--- a/input/weight
+++ b/input/weight
@@ -1,37 +1,76 @@
-1  10 
-S   8   1.00
-  24350.0000000              0.0005020        
-   3650.0000000              0.0038810        
-    829.6000000              0.0199970        
-    234.0000000              0.0784180        
-     75.6100000              0.2296760        
-     26.7300000              0.4327220        
-      9.9270000              0.3506420        
-      1.1020000             -0.0076450        
-S   8   1.00
-  24350.0000000             -0.0001180        
-   3650.0000000             -0.0009150        
-    829.6000000             -0.0047370        
-    234.0000000             -0.0192330        
-     75.6100000             -0.0603690        
-     26.7300000             -0.1425080        
-      9.9270000             -0.1777100        
-      1.1020000              0.6058360        
+1     9
+S   9   1.00
+      1.471000D+04           7.210000D-04
+      2.207000D+03           5.553000D-03
+      5.028000D+02           2.826700D-02
+      1.426000D+02           1.064440D-01
+      4.647000D+01           2.868140D-01
+      1.670000D+01           4.486410D-01
+      6.356000D+00           2.647610D-01
+      1.316000D+00           1.533300D-02
+      3.897000D-01          -2.332000D-03
+S   9   1.00
+      1.471000D+04          -1.650000D-04
+      2.207000D+03          -1.308000D-03
+      5.028000D+02          -6.495000D-03
+      1.426000D+02          -2.669100D-02
+      4.647000D+01          -7.369000D-02
+      1.670000D+01          -1.707760D-01
+      6.356000D+00          -1.123270D-01
+      1.316000D+00           5.628140D-01
+      3.897000D-01           5.687780D-01
 S   1   1.00
-      2.8360000              1.0000000        
+      3.897000D-01           1.000000D+00
 S   1   1.00
-      0.3782000              1.0000000        
-P   3   1.00
-     54.7000000              0.0171510        
-     12.4300000              0.1076560        
-      3.6790000              0.3216810        
+      0.0986300              1.0000000
+P   4   1.00
+      2.267000D+01           4.487800D-02
+      4.977000D+00           2.357180D-01
+      1.347000D+00           5.085210D-01
+      3.471000D-01           4.581200D-01
 P   1   1.00
-      1.1430000              1.0000000        
+      3.471000D-01           1.000000D+00
 P   1   1.00
-      0.3300000              1.0000000        
+      0.0850200              1.0000000
 D   1   1.00
-      4.0140000              1.0000000        
+      1.640000D+00           1.0000000
 D   1   1.00
-      1.0960000              1.0000000        
-F   1   1.00
-      2.5440000              1.0000000   
+      0.4640000              1.0000000
+2   9
+S   9   1.00
+      1.471000D+04           7.210000D-04
+      2.207000D+03           5.553000D-03
+      5.028000D+02           2.826700D-02
+      1.426000D+02           1.064440D-01
+      4.647000D+01           2.868140D-01
+      1.670000D+01           4.486410D-01
+      6.356000D+00           2.647610D-01
+      1.316000D+00           1.533300D-02
+      3.897000D-01          -2.332000D-03
+S   9   1.00
+      1.471000D+04          -1.650000D-04
+      2.207000D+03          -1.308000D-03
+      5.028000D+02          -6.495000D-03
+      1.426000D+02          -2.669100D-02
+      4.647000D+01          -7.369000D-02
+      1.670000D+01          -1.707760D-01
+      6.356000D+00          -1.123270D-01
+      1.316000D+00           5.628140D-01
+      3.897000D-01           5.687780D-01
+S   1   1.00
+      3.897000D-01           1.000000D+00
+S   1   1.00
+      0.0986300              1.0000000
+P   4   1.00
+      2.267000D+01           4.487800D-02
+      4.977000D+00           2.357180D-01
+      1.347000D+00           5.085210D-01
+      3.471000D-01           4.581200D-01
+P   1   1.00
+      3.471000D-01           1.000000D+00
+P   1   1.00
+      0.0850200              1.0000000
+D   1   1.00
+      1.640000D+00           1.0000000
+D   1   1.00
+      0.4640000              1.0000000
diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90
index 6f8b47e..301b31d 100644
--- a/src/QuAcK/G0T0.f90
+++ b/src/QuAcK/G0T0.f90
@@ -1,4 +1,4 @@
-subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,nVV,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0T0)
+subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF,eG0T0)
 
 ! Perform G0W0 calculation with a T-matrix self-energy (G0T0)
 
@@ -12,25 +12,19 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,n
   logical,intent(in)            :: triplet_manifold
   double precision,intent(in)   :: eta
 
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nOO,nVV
+  integer,intent(in)            :: nBas,nC,nO,nV,nR
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: ERHF
   double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: cHF(nBas,nBas)
-  double precision,intent(in)   :: PHF(nBas,nBas)
-  double precision,intent(in)   :: Hc(nBas,nBas)
-  double precision,intent(in)   :: H(nBas,nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
 ! Local variables
 
-  logical                       :: dRPA
-  integer                       :: ispin,jspin
+  integer                       :: ispin
+  integer                       :: nOO
+  integer                       :: nVV
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
-  double precision              :: EcGM
-  double precision,allocatable  :: SigT(:)
-  double precision,allocatable  :: Z(:)
   double precision,allocatable  :: Omega1(:,:)
   double precision,allocatable  :: X1(:,:,:)
   double precision,allocatable  :: Y1(:,:,:)
@@ -39,6 +33,8 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,n
   double precision,allocatable  :: X2(:,:,:)
   double precision,allocatable  :: Y2(:,:,:)
   double precision,allocatable  :: rho2(:,:,:,:)
+  double precision,allocatable  :: SigT(:)
+  double precision,allocatable  :: Z(:)
 
 ! Output variables
 
@@ -56,6 +52,9 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,n
 
   ispin = 1
 
+  nOO = nO*(nO+1)/2
+  nVV = nV*(nV+1)/2
+
 ! Memory allocation
 
   allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & 
@@ -65,29 +64,29 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,n
 
 ! Compute linear response
 
-  call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,eHF,ERI, & 
-                          Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),    & 
-                          Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),    & 
+  call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,eHF(:),ERI(:,:,:,:), & 
+                          Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),                & 
+                          Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),                & 
                           EcRPA(ispin))
 
 ! Compute excitation densities for the T-matrix
 
-  call excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,                     & 
+  call excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI(:,:,:,:),            & 
                                   X1(:,:,ispin),Y1(:,:,ispin),rho1(:,:,:,ispin), & 
                                   X2(:,:,ispin),Y2(:,:,ispin),rho2(:,:,:,ispin))
 
 ! Compute T-matrix version of the self-energy 
 
-  call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,  & 
-                                Omega1(:,ispin),rho1(:,:,:,ispin), &  
-                                Omega2(:,ispin),rho2(:,:,:,ispin), & 
-                                SigT)
+  call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF(:), & 
+                                Omega1(:,ispin),rho1(:,:,:,ispin),   &  
+                                Omega2(:,ispin),rho2(:,:,:,ispin),   & 
+                                SigT(:))
 
 ! Compute renormalization factor for T-matrix self-energy
 
-  call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,  & 
-                                      Omega1(:,ispin),rho1(:,:,:,ispin), & 
-                                      Omega2(:,ispin),rho2(:,:,:,ispin), & 
+  call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF(:), & 
+                                      Omega1(:,ispin),rho1(:,:,:,ispin),   & 
+                                      Omega2(:,ispin),rho2(:,:,:,ispin),   & 
                                       Z(:))
 
 ! Solve the quasi-particle equation
@@ -99,7 +98,7 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,n
   call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
   call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
 
-  call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcRPA(ispin))
+  call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(ispin))
 
 ! Perform BSE calculation
 
diff --git a/src/QuAcK/GF2_diag.f90 b/src/QuAcK/GF2_diag.f90
index 18e323e..28ddb63 100644
--- a/src/QuAcK/GF2_diag.f90
+++ b/src/QuAcK/GF2_diag.f90
@@ -1,4 +1,4 @@
-subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
+subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
 ! Perform second-order Green function calculation in diagonal approximation
 
@@ -10,6 +10,7 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
   integer,intent(in)            :: maxSCF
   double precision,intent(in)   :: thresh
   integer,intent(in)            :: max_diis
+  logical,intent(in)            :: linearize
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nO
   integer,intent(in)            :: nC
@@ -28,6 +29,7 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
   double precision,allocatable  :: eGF2(:)
   double precision,allocatable  :: eOld(:)
   double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Z(:)
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: e_diis(:,:)
 
@@ -43,7 +45,7 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
 
 ! Memory allocation
 
-  allocate(Bpp(nBas,2),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
+  allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
 
 ! Initialization
 
@@ -95,7 +97,49 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
       end do
     end do
 
-    eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
+    ! 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
 
     Conv = maxval(abs(eGF2 - eOld))
 
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index af170fd..928c5b0 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -16,7 +16,7 @@ program QuAcK
 
   integer                       :: nNuc,nBas,nBasCABS
   integer                       :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
-  integer                       :: nS(nspin),nOO(nspin),nVV(nspin)
+  integer                       :: nS(nspin)
   double precision              :: ENuc,ERHF,EUHF,Norm
   double precision              :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
 
@@ -446,7 +446,7 @@ program QuAcK
   if(doGF2) then
 
     call cpu_time(start_GF2)
-    call GF2_diag(maxSCF_GF,thresh_GF,n_diis_GF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
+    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 cpu_time(end_GF2)
 
     t_GF2 = end_GF2 - start_GF2
@@ -533,12 +533,9 @@ program QuAcK
 
   if(doG0T0) then
     
-    nOO(:) = nO(:)*(nO(:)+1)/2
-    nVV(:) = nV(:)*(nV(:)+1)/2
-
     call cpu_time(start_G0T0)
     call G0T0(BSE,singlet_manifold,triplet_manifold,eta, & 
-              nBas,nC(1),nO(1),nV(1),nR(1),nOO(1),nVV(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0T0)
+              nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
     call cpu_time(end_G0T0)
   
     t_G0T0 = end_G0T0 - start_G0T0
diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90
index 530623a..9b594d6 100644
--- a/src/QuAcK/excitation_density_Tmatrix.f90
+++ b/src/QuAcK/excitation_density_Tmatrix.f90
@@ -39,7 +39,8 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2
          do d=nO+1,c
             cd = cd + 1
             rho1(p,i,ab) = rho1(p,i,ab) & 
-                         + (ERI(p,i,c,d) - 0.5d0*ERI(p,i,d,c))*X1(cd,ab)!/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
+                         + 1.0d0*(ERI(p,i,c,d) - 1.0d0*ERI(p,i,d,c))*X1(cd,ab)!& 
+!                        /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
           enddo
         enddo
 
@@ -48,7 +49,8 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2
          do l=nC+1,k
             kl = kl + 1
             rho1(p,i,ab) = rho1(p,i,ab) & 
-                         + (ERI(p,i,k,l) - 0.5d0*ERI(p,i,l,k))*Y1(kl,ab)!/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
+                         + 1.0d0*(ERI(p,i,k,l) - 1.0d0*ERI(p,i,l,k))*Y1(kl,ab)!& 
+!                        /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
           enddo
         enddo
 
@@ -63,7 +65,8 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2
          do d=nO+1,c
             cd = cd + 1
             rho2(p,a,ij) = rho2(p,a,ij) & 
-                         + (ERI(p,a,c,d) - 0.5d0*ERI(p,a,d,c))*X2(cd,ij)!/sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d)))
+                         + 1.0d0*(ERI(p,a,c,d) - 1.0d0*ERI(p,a,d,c))*X2(cd,ij)!&
+!                        /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d)))
           enddo
         enddo
 
@@ -72,7 +75,8 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2
          do l=nC+1,k
             kl = kl + 1
             rho2(p,a,ij) = rho2(p,a,ij) & 
-                         + (ERI(p,a,k,l) - 0.5d0*ERI(p,a,l,k))*Y2(kl,ij)!/sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l)))
+                         + 1.0d0*(ERI(p,a,k,l) - 1.0d0*ERI(p,a,l,k))*Y2(kl,ij)!& 
+!                        /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l)))
           enddo
         enddo
 
diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90
index 6671e15..c566837 100644
--- a/src/QuAcK/linear_response_pp.f90
+++ b/src/QuAcK/linear_response_pp.f90
@@ -60,8 +60,8 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
 
 ! Off-diagonal blocks
 
-  M(    1:nVV    ,nVV+1:nOO+nVV) = -           B(1:nVV,1:nOO)
-  M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(B(1:nVV,1:nOO))
+  M(    1:nVV    ,nVV+1:nOO+nVV) = +           B(1:nVV,1:nOO)
+  M(nVV+1:nOO+nVV,    1:nVV)     = - transpose(B(1:nVV,1:nOO))
 
 ! print*, 'pp-RPA matrix'
 ! call matout(nOO+nVV,nOO+nVV,M(:,:))
@@ -75,19 +75,9 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
 ! call matout(nOO+nVV,1,Omega(:))
 ! write(*,*) 
 
-! call matout(nOO+nVV,nOO+nVV,matmul(transpose(Z(:,:)),matmul(W(:,:),Z(:,:))))
-! write(*,*) 
-
 ! Split the various quantities in p-p and h-h parts
 
   call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
-! Omega1(:) = Omega(nOO+1:nOO+nVV)
-! Omega2(:) = Omega(1:nOO)
-
-! X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV)
-! Y1(:,:) = Z(    1:nOO    ,nOO+1:nOO+nVV)
-! X2(:,:) = Z(nOO+1:nOO+nVV,    1:nOO    ) 
-! Y2(:,:) = Z(    1:nOO    ,nOO+1:nOO+nVV)
 
 ! Compute the RPA correlation energy
 
@@ -95,6 +85,15 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
   Ec_ppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
 ! Ec_ppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
 
-  print*,'Ec(pp-RPA) = ',Ec_ppRPA
+! write(*,*)'X1'
+! call matout(nVV,nVV,X1)
+! write(*,*)'Y1'
+! call matout(nVV,nOO,Y1)
+! write(*,*)'X2'
+! call matout(nOO,nVV,X2)
+! write(*,*)'Y2'
+! call matout(nOO,nOO,Y2)
+
+! print*,'Ec(pp-RPA) = ',Ec_ppRPA
 
 end subroutine linear_response_pp
diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90
index 112982b..00681ba 100644
--- a/src/QuAcK/sort_ppRPA.f90
+++ b/src/QuAcK/sort_ppRPA.f90
@@ -15,7 +15,6 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 ! Local variables
 
   integer                       :: pq,ab,ij
-  double precision,allocatable  :: s(:,:)
 
 ! Output variables
 
@@ -26,10 +25,6 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
   double precision,intent(out)  :: X2(nVV,nOO)
   double precision,intent(out)  :: Y2(nOO,nOO)
 
-! Memory allocation
-
-  allocate(s(nOO+nVV,nOO+nVV))
-
 ! Initializatiom
 
   Omega1(:) = 0d0
@@ -68,21 +63,4 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 ! call matout(nOO,1,Omega2(:))
 ! write(*,*)
 
-! Check eigenvector signatures
-
-! s(    1:    nVV,    1:    nVV) = matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1)
-! s(nVV+1:nOO+nVV,nVV+1:nOO+nVV) = matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2)
-
-! write(*,*) 'Signatures pp'
-! do ab=1,nVV
-!   write(*,'(I6,F10.6)') ab,s(ab,ab)
-! end do
-! write(*,*)
-
-! write(*,*) 'Signatures hh'
-! do ij=1,nOO
-!   write(*,'(I6,F10.6)') ij,s(ij,ij)
-! end do
-! write(*,*)
-
 end subroutine sort_ppRPA