diff --git a/input/basis b/input/basis
index b18e129..35be48a 100644
--- a/input/basis
+++ b/input/basis
@@ -1,52 +1,39 @@
-1     9
+1  10
 S   9   1.00
-      1.469000D+03           7.660000D-04
-      2.205000D+02           5.892000D-03
-      5.026000D+01           2.967100D-02
-      1.424000D+01           1.091800D-01
-      4.581000D+00           2.827890D-01
-      1.580000D+00           4.531230D-01
-      5.640000D-01           2.747740D-01
-      7.345000D-02           9.751000D-03
-      2.805000D-02          -3.180000D-03
+   6863.0000000              0.0002360        
+   1030.0000000              0.0018260        
+    234.7000000              0.0094520        
+     66.5600000              0.0379570        
+     21.6900000              0.1199650        
+      7.7340000              0.2821620        
+      2.9160000              0.4274040        
+      1.1300000              0.2662780        
+      0.1101000             -0.0072750        
 S   9   1.00
-      1.469000D+03          -1.200000D-04
-      2.205000D+02          -9.230000D-04
-      5.026000D+01          -4.689000D-03
-      1.424000D+01          -1.768200D-02
-      4.581000D+00          -4.890200D-02
-      1.580000D+00          -9.600900D-02
-      5.640000D-01          -1.363800D-01
-      7.345000D-02           5.751020D-01
-      2.805000D-02           5.176610D-01
+   6863.0000000             -0.0000430        
+   1030.0000000             -0.0003330        
+    234.7000000             -0.0017360        
+     66.5600000             -0.0070120        
+     21.6900000             -0.0231260        
+      7.7340000             -0.0581380        
+      2.9160000             -0.1145560        
+      1.1300000             -0.1359080        
+      0.1101000              0.5774410        
 S   1   1.00
-      2.805000D-02           1.000000D+00
+      0.2577000              1.0000000        
 S   1   1.00
-      8.600000D-03           1.000000D+00
-P   4   1.00
-      1.534000D+00           2.278400D-02
-      2.749000D-01           1.391070D-01
-      7.362000D-02           5.003750D-01
-      2.403000D-02           5.084740D-01
+      0.0440900              1.0000000        
+P   3   1.00
+      7.4360000              0.0107360        
+      1.5770000              0.0628540        
+      0.4352000              0.2481800        
 P   1   1.00
-      2.403000D-02           1.000000D+00
+      0.1438000              1.0000000        
 P   1   1.00
-      5.800000D-03           1.000000D+00
+      0.0499400              1.0000000        
 D   1   1.00
-      1.144000D-01           1.0000000
+      0.3480000              1.0000000        
 D   1   1.00
-      7.330000D-02           1.000000D+00
-2     5
-S   4   1.00
-      1.301000D+01           1.968500D-02
-      1.962000D+00           1.379770D-01
-      4.446000D-01           4.781480D-01
-      1.220000D-01           5.012400D-01
-S   1   1.00
-      1.220000D-01           1.000000D+00
-S   1   1.00
-      0.0297400              1.0000000
-P   1   1.00
-      7.270000D-01           1.0000000
-P   1   1.00
-      0.1410000              1.0000000
+      0.1803000              1.0000000        
+F   1   1.00
+      0.3250000              1.0000000   
diff --git a/input/methods b/input/methods
index 513d614..9ac4f33 100644
--- a/input/methods
+++ b/input/methods
@@ -7,7 +7,7 @@
 # CIS TDHF ppRPA ADC
   F   F    F     F
 # GF2 GF3
-  T   F
+  F   F
 # G0W0 evGW qsGW 
   F    F    F
 # G0T0 evGT qsGT 
diff --git a/input/molecule b/input/molecule
index c38d7a8..6a6f6d1 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,5 +1,4 @@
 # nAt nEla nElb nCore nRyd
- 2   2     2   0     0
+  1   2     2   0     0
 # Znuc   x            y           z
-  Li    0.           0.         0.
-  H     0.           0.         3.09839495
+  Be      0.0          0.0         0.0
diff --git a/input/weight b/input/weight
index b18e129..35be48a 100644
--- a/input/weight
+++ b/input/weight
@@ -1,52 +1,39 @@
-1     9
+1  10
 S   9   1.00
-      1.469000D+03           7.660000D-04
-      2.205000D+02           5.892000D-03
-      5.026000D+01           2.967100D-02
-      1.424000D+01           1.091800D-01
-      4.581000D+00           2.827890D-01
-      1.580000D+00           4.531230D-01
-      5.640000D-01           2.747740D-01
-      7.345000D-02           9.751000D-03
-      2.805000D-02          -3.180000D-03
+   6863.0000000              0.0002360        
+   1030.0000000              0.0018260        
+    234.7000000              0.0094520        
+     66.5600000              0.0379570        
+     21.6900000              0.1199650        
+      7.7340000              0.2821620        
+      2.9160000              0.4274040        
+      1.1300000              0.2662780        
+      0.1101000             -0.0072750        
 S   9   1.00
-      1.469000D+03          -1.200000D-04
-      2.205000D+02          -9.230000D-04
-      5.026000D+01          -4.689000D-03
-      1.424000D+01          -1.768200D-02
-      4.581000D+00          -4.890200D-02
-      1.580000D+00          -9.600900D-02
-      5.640000D-01          -1.363800D-01
-      7.345000D-02           5.751020D-01
-      2.805000D-02           5.176610D-01
+   6863.0000000             -0.0000430        
+   1030.0000000             -0.0003330        
+    234.7000000             -0.0017360        
+     66.5600000             -0.0070120        
+     21.6900000             -0.0231260        
+      7.7340000             -0.0581380        
+      2.9160000             -0.1145560        
+      1.1300000             -0.1359080        
+      0.1101000              0.5774410        
 S   1   1.00
-      2.805000D-02           1.000000D+00
+      0.2577000              1.0000000        
 S   1   1.00
-      8.600000D-03           1.000000D+00
-P   4   1.00
-      1.534000D+00           2.278400D-02
-      2.749000D-01           1.391070D-01
-      7.362000D-02           5.003750D-01
-      2.403000D-02           5.084740D-01
+      0.0440900              1.0000000        
+P   3   1.00
+      7.4360000              0.0107360        
+      1.5770000              0.0628540        
+      0.4352000              0.2481800        
 P   1   1.00
-      2.403000D-02           1.000000D+00
+      0.1438000              1.0000000        
 P   1   1.00
-      5.800000D-03           1.000000D+00
+      0.0499400              1.0000000        
 D   1   1.00
-      1.144000D-01           1.0000000
+      0.3480000              1.0000000        
 D   1   1.00
-      7.330000D-02           1.000000D+00
-2     5
-S   4   1.00
-      1.301000D+01           1.968500D-02
-      1.962000D+00           1.379770D-01
-      4.446000D-01           4.781480D-01
-      1.220000D-01           5.012400D-01
-S   1   1.00
-      1.220000D-01           1.000000D+00
-S   1   1.00
-      0.0297400              1.0000000
-P   1   1.00
-      7.270000D-01           1.0000000
-P   1   1.00
-      0.1410000              1.0000000
+      0.1803000              1.0000000        
+F   1   1.00
+      0.3250000              1.0000000   
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 928c5b0..b346369 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -536,6 +536,7 @@ program QuAcK
     call cpu_time(start_G0T0)
     call G0T0(BSE,singlet_manifold,triplet_manifold,eta, & 
               nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
+    call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
     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 751690d..65a6d66 100644
--- a/src/QuAcK/excitation_density_Tmatrix.f90
+++ b/src/QuAcK/excitation_density_Tmatrix.f90
@@ -46,9 +46,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
            do d=nO+1,c
               cd = cd + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
-                           + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) 
-!                          + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
-!                          /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
+!                          + ERI(p,i,c,d)*X1(cd,ab) 
+                           + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
+                           /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
             end do
           end do
 
@@ -57,9 +57,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
            do l=nC+1,k
               kl = kl + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
-                           + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
-!                          + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
-!                          /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
+!                          + ERI(p,i,k,l)*Y1(kl,ab)
+                           + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
+                           /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
             end do
           end do
 
@@ -74,9 +74,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
            do d=nO+1,c
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij)
-!                          + (ERI(p,a,c,d) + ERI(p,a,d,c))*X2(cd,ij) &
-!                          /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d)))
+!                          + ERI(p,a,c,d)*X2(cd,ij)
+                           + (ERI(p,a,c,d) + ERI(p,a,d,c))*X2(cd,ij) &
+                           /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d)))
 
             end do
           end do
@@ -86,9 +86,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
            do l=nC+1,k
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)  
-!                          + (ERI(p,a,k,l) + ERI(p,a,l,k))*Y2(kl,ij) &
-!                          /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l)))
+!                          + ERI(p,a,k,l)*Y2(kl,ij)  
+                           + (ERI(p,a,k,l) + ERI(p,a,l,k))*Y2(kl,ij) &
+                           /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l)))
 
             end do
           end do
diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90
index c566837..3745275 100644
--- a/src/QuAcK/linear_response_pp.f90
+++ b/src/QuAcK/linear_response_pp.f90
@@ -1,4 +1,4 @@
-subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA)
+subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
 
 ! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013)
 
@@ -32,7 +32,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
   double precision,intent(out)  :: Omega2(nOO)
   double precision,intent(out)  :: X2(nVV,nOO)
   double precision,intent(out)  :: Y2(nOO,nOO)
-  double precision,intent(out)  :: Ec_ppRPA
+  double precision,intent(out)  :: EcppRPA
 
 ! Memory allocation
 
@@ -81,9 +81,9 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
 
 ! Compute the RPA correlation energy
 
-! Ec_ppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
-  Ec_ppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
-! Ec_ppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
+! EcppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
+  EcppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
+! EcppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
 
 ! write(*,*)'X1'
 ! call matout(nVV,nVV,X1)
@@ -94,6 +94,6 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
 ! write(*,*)'Y2'
 ! call matout(nOO,nOO,Y2)
 
-! print*,'Ec(pp-RPA) = ',Ec_ppRPA
+! print*,'Ec(pp-RPA) = ',EcppRPA
 
 end subroutine linear_response_pp
diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90
index 118ac64..76a0c5a 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix.f90
@@ -11,8 +11,8 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
 
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas,nC,nO,nV,nR
-  integer,intent(in)            :: nOOs,nVVs
-  integer,intent(in)            :: nOOt,nVVt
+  integer,intent(in)            :: nOOs,nOOt
+  integer,intent(in)            :: nVVs,nVVt
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Omega1s(nVVs),Omega1t(nVVt)
   double precision,intent(in)   :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt)
@@ -100,40 +100,6 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     enddo
   enddo
 
-!----------------------------------------------
-! Singlet part of the T-matrix self-energy
-!----------------------------------------------
-
-! Occupied part of the T-matrix self-energy 
-
-  do p=nC+1,nBas-nR
-    do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=nO+1,c
-          cd = cd + 1
-          eps = e(p) + e(i) - Omega1s(cd)
-          Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2
-        enddo
-      enddo
-    enddo
-  enddo
-
-! Virtual part of the T-matrix self-energy
-
-  do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=nC+1,k
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2s(kl)
-          Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2
-        enddo
-      enddo
-    enddo
-  enddo
-
 ! Compute renormalization factor from derivative of SigT
  
   Z(:) = 1d0/(1d0 - Z(:))
diff --git a/src/QuAcK/renormalization_factor_Tmatrix_so.f90 b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
new file mode 100644
index 0000000..e9eb7dd
--- /dev/null
+++ b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
@@ -0,0 +1,71 @@
+subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z)
+
+! Compute renormalization factor of the T-matrix self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas,nC,nO,nV,nR
+  integer,intent(in)            :: nOO
+  integer,intent(in)            :: nVV
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Omega1(nVV)
+  double precision,intent(in)   :: rho1(nBas,nBas,nVV)
+  double precision,intent(in)   :: Omega2(nOO)
+  double precision,intent(in)   :: rho2(nBas,nBas,nOO)
+
+! Local variables
+
+  integer                       :: i,j,k,l,a,b,c,d,p,cd,kl
+  double precision              :: eps
+
+! Output variables
+
+  double precision,intent(out)  :: Z(nBas)
+
+! Initialize
+
+  Z(:)  = 0d0
+
+!----------------------------------------------
+! T-matrix renormalization factor in the spinorbital basis
+!----------------------------------------------
+
+! Occupied part of the T-matrix self-energy 
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      cd = 0
+      do c=nO+1,nBas-nR
+        do d=nO+1,c-1
+          cd = cd + 1
+          eps = e(p) + e(i) - Omega1(cd)
+          Z(p) = Z(p) - rho1(p,i,cd)**2/eps**2
+        enddo
+      enddo
+    enddo
+  enddo
+
+! Virtual part of the T-matrix self-energy
+
+  do p=nC+1,nBas-nR
+    do a=nO+1,nBas-nR
+      kl = 0
+      do k=nC+1,nO
+        do l=nC+1,k-1
+          kl = kl + 1
+          eps = e(p) + e(a) - Omega2(kl)
+          Z(p) = Z(p) - rho2(p,a,kl)**2/eps**2
+        enddo
+      enddo
+    enddo
+  enddo
+
+! Compute renormalization factor from derivative of SigT
+ 
+  Z(:) = 1d0/(1d0 - Z(:))
+
+end subroutine renormalization_factor_Tmatrix_so
diff --git a/src/QuAcK/self_energy_Tmatrix_diag_so.f90 b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
new file mode 100644
index 0000000..c35dedd
--- /dev/null
+++ b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
@@ -0,0 +1,71 @@
+subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT)
+
+! Compute diagonal of the correlation part of the T-matrix self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nOO
+  integer,intent(in)            :: nVV
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Omega1(nVV)
+  double precision,intent(in)   :: rho1(nBas,nBas,nVV)
+  double precision,intent(in)   :: Omega2(nOO)
+  double precision,intent(in)   :: rho2(nBas,nBas,nOO)
+
+! Local variables
+
+  integer                       :: i,j,k,l,a,b,c,d,p,cd,kl
+  double precision              :: eps
+
+! Output variables
+
+  double precision,intent(out)  :: SigT(nBas)
+
+! Initialize 
+
+  SigT(:) = 0d0
+
+!----------------------------------------------
+! T-matrix self-energy in the spinorbital basis
+!----------------------------------------------
+
+! Occupied part of the T-matrix self-energy 
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      cd = 0
+      do c=nO+1,nBas-nR
+        do d=nO+1,c-1
+          cd = cd + 1
+          eps = e(p) + e(i) - Omega1(cd)
+          SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
+        enddo
+      enddo
+    enddo
+  enddo
+
+! Virtual part of the T-matrix self-energy
+
+  do p=nC+1,nBas-nR
+    do a=nO+1,nBas-nR
+      kl = 0
+      do k=nC+1,nO
+        do l=nC+1,k-1
+          kl = kl + 1
+          eps = e(p) + e(a) - Omega2(kl)
+          SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine self_energy_Tmatrix_diag_so
diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90
new file mode 100644
index 0000000..5543682
--- /dev/null
+++ b/src/QuAcK/soG0T0.f90
@@ -0,0 +1,127 @@
+subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
+
+! Perform G0W0 calculation with a T-matrix self-energy (G0T0) in the spinorbital basis
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+
+  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)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: ispin
+  integer                       :: nOO
+  integer                       :: nVV
+  double precision              :: EcRPA
+  integer                       :: nBas2,nC2,nO2,nV2,nR2
+  double precision,allocatable  :: Omega1(:)
+  double precision,allocatable  :: X1(:,:)
+  double precision,allocatable  :: Y1(:,:)
+  double precision,allocatable  :: rho1(:,:,:)
+  double precision,allocatable  :: Omega2(:)
+  double precision,allocatable  :: X2(:,:)
+  double precision,allocatable  :: Y2(:,:)
+  double precision,allocatable  :: rho2(:,:,:)
+  double precision,allocatable  :: SigT(:)
+  double precision,allocatable  :: Z(:)
+  double precision,allocatable  :: eG0T0(:)
+  double precision,allocatable  :: seHF(:)
+  double precision,allocatable  :: sERI(:,:,:,:)
+
+! Output variables
+
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|          One-shot soG0T0 calculation         |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Spatial to spin orbitals
+
+  nBas2 = 2*nBas
+
+  allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
+
+  call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
+  call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
+
+! Define occupied and virtual spaces
+
+  nO2 = 2*nO
+  nV2 = 2*nV
+  nC2 = 2*nC 
+  nR2 = 2*nR
+
+! Dimensions of the rr-RPA linear reponse matrices
+
+  nOO = nO2*(nO2-1)/2
+  nVV = nV2*(nV2-1)/2
+
+! Memory allocation
+
+  allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & 
+           Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & 
+           rho1(nBas2,nBas2,nVV),rho2(nBas2,nBas2,nOO), & 
+           eG0T0(nBas2),SigT(nBas2),Z(nBas2))
+
+!----------------------------------------------
+! Spinorbital basis
+!----------------------------------------------
+
+ ispin = 2
+
+! Compute linear response
+
+  call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2, & 
+                          nOO,nVV,seHF(:),sERI(:,:,:,:),  & 
+                          Omega1(:),X1(:,:),Y1(:,:),   & 
+                          Omega2(:),X2(:,:),Y2(:,:),   & 
+                          EcRPA)
+
+  call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
+  call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:))
+
+! Compute excitation densities for the T-matrix
+
+  call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nR2,nOO,nVV,sERI(:,:,:,:), & 
+                                  X1(:,:),Y1(:,:),rho1(:,:,:),             & 
+                                  X2(:,:),Y2(:,:),rho2(:,:,:))
+
+!----------------------------------------------
+! Compute T-matrix version of the self-energy 
+!----------------------------------------------
+
+  call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
+                                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(:,:,:), & 
+                                      Z(:))
+
+!----------------------------------------------
+! Solve the quasi-particle equation
+!----------------------------------------------
+
+  eG0T0(:) = seHF(:) + Z(:)*SigT(:)
+
+!----------------------------------------------
+! Dump results
+!----------------------------------------------
+
+  call print_G0T0(nBas2,nO2,seHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA)
+
+end subroutine soG0T0
diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90
index 00681ba..11951e8 100644
--- a/src/QuAcK/sort_ppRPA.f90
+++ b/src/QuAcK/sort_ppRPA.f90
@@ -32,6 +32,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
   ab = 0
   ij = 0
+
   do pq=1,nOO+nVV
 
     if(Omega(pq) > 0d0) then 
diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90
index 700aee4..61edacf 100644
--- a/src/utils/wrap_lapack.f90
+++ b/src/utils/wrap_lapack.f90
@@ -53,7 +53,7 @@ subroutine diagonalize_general_matrix(N,A,e,X)
   lwork = 4*N
   allocate(WI(N),VL(N,N),work(lwork))
 
-  call dgeev('V','V',N,A,N,e,WI,VL,N,X,N,work,lwork,info)
+  call dgeev('N','V',N,A,N,e,WI,VL,N,X,N,work,lwork,info)
  
   if(info /= 0) then 
     print*,'Problem in diagonalize_matrix (dgeev)!!'