diff --git a/input/basis b/input/basis
index 35be48a..b18e129 100644
--- a/input/basis
+++ b/input/basis
@@ -1,39 +1,52 @@
-1  10
+1     9
 S   9   1.00
-   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        
+      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
 S   9   1.00
-   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        
+      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
 S   1   1.00
-      0.2577000              1.0000000        
+      2.805000D-02           1.000000D+00
 S   1   1.00
-      0.0440900              1.0000000        
-P   3   1.00
-      7.4360000              0.0107360        
-      1.5770000              0.0628540        
-      0.4352000              0.2481800        
+      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
 P   1   1.00
-      0.1438000              1.0000000        
+      2.403000D-02           1.000000D+00
 P   1   1.00
-      0.0499400              1.0000000        
+      5.800000D-03           1.000000D+00
 D   1   1.00
-      0.3480000              1.0000000        
+      1.144000D-01           1.0000000
 D   1   1.00
-      0.1803000              1.0000000        
-F   1   1.00
-      0.3250000              1.0000000   
+      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
diff --git a/input/methods b/input/methods
index 9ac4f33..6264e2f 100644
--- a/input/methods
+++ b/input/methods
@@ -5,7 +5,7 @@
 # CCD CCSD CCSD(T)
   F   F    F
 # CIS TDHF ppRPA ADC
-  F   F    F     F
+  F   F    T     F
 # GF2 GF3
   F   F
 # G0W0 evGW qsGW 
diff --git a/input/molecule b/input/molecule
index 6a6f6d1..c38d7a8 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,5 @@
 # nAt nEla nElb nCore nRyd
-  1   2     2   0     0
+ 2   2     2   0     0
 # Znuc   x            y           z
-  Be      0.0          0.0         0.0
+  Li    0.           0.         0.
+  H     0.           0.         3.09839495
diff --git a/input/weight b/input/weight
index 35be48a..b18e129 100644
--- a/input/weight
+++ b/input/weight
@@ -1,39 +1,52 @@
-1  10
+1     9
 S   9   1.00
-   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        
+      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
 S   9   1.00
-   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        
+      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
 S   1   1.00
-      0.2577000              1.0000000        
+      2.805000D-02           1.000000D+00
 S   1   1.00
-      0.0440900              1.0000000        
-P   3   1.00
-      7.4360000              0.0107360        
-      1.5770000              0.0628540        
-      0.4352000              0.2481800        
+      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
 P   1   1.00
-      0.1438000              1.0000000        
+      2.403000D-02           1.000000D+00
 P   1   1.00
-      0.0499400              1.0000000        
+      5.800000D-03           1.000000D+00
 D   1   1.00
-      0.3480000              1.0000000        
+      1.144000D-01           1.0000000
 D   1   1.00
-      0.1803000              1.0000000        
-F   1   1.00
-      0.3250000              1.0000000   
+      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
diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90
index 65a6d66..a1506dd 100644
--- a/src/QuAcK/excitation_density_Tmatrix.f90
+++ b/src/QuAcK/excitation_density_Tmatrix.f90
@@ -158,4 +158,63 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
 
   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=nO+1,nBas-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,a,c,d) - ERI(p,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,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)
+            end do
+          end do
+
+        end do
+      end do
+    end do
+
+  end if
+
 end subroutine excitation_density_Tmatrix
diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90
index 696b490..e9d18a5 100644
--- a/src/QuAcK/linear_response_B_pp.f90
+++ b/src/QuAcK/linear_response_B_pp.f90
@@ -64,4 +64,26 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 
   end if
 
+! Build the B matrix for the spinorbital basis
+
+  if(ispin == 3) then
+
+    ab = 0
+    do a=nO+1,nBas-nR
+     do b=a+1,nBas-nR
+        ab = ab + 1
+        ij = 0
+        do i=nC+1,nO
+         do j=i+1,nO
+            ij = ij + 1
+ 
+            B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
+
 end subroutine linear_response_B_pp
diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90
index c68b257..72e56f2 100644
--- a/src/QuAcK/linear_response_C_pp.f90
+++ b/src/QuAcK/linear_response_C_pp.f90
@@ -72,4 +72,27 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 
   end if
 
+! Build C matrix for the spinorbital basis
+
+  if(ispin == 3) then
+
+    ab = 0
+    do a=nO+1,nBas-nR
+     do b=a+1,nBas-nR
+        ab = ab + 1
+        cd = 0
+        do c=nO+1,nBas-nR
+         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) & 
+                          + ERI(a,b,c,d) - ERI(a,b,d,c)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
+
 end subroutine linear_response_C_pp
diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90
index a686d6e..fec3d10 100644
--- a/src/QuAcK/linear_response_D_pp.f90
+++ b/src/QuAcK/linear_response_D_pp.f90
@@ -72,4 +72,27 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 
   end if
 
+! Build the D matrix for the spinorbital basis
+
+  if(ispin == 3) then
+
+    ij = 0
+    do i=nC+1,nO
+     do j=i+1,nO
+        ij = ij + 1
+        kl = 0
+        do k=nC+1,nO
+         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) & 
+                          + ERI(i,j,k,l) - ERI(i,j,l,k)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
+
 end subroutine linear_response_D_pp
diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90
index d8e8a60..fca9912 100644
--- a/src/QuAcK/print_excitation.f90
+++ b/src/QuAcK/print_excitation.f90
@@ -19,6 +19,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
 
   if(ispin == 1) spin_manifold = 'singlet'
   if(ispin == 2) spin_manifold = 'triplet'
+  if(ispin == 3) spin_manifold = 'spinorb'
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------'
diff --git a/src/QuAcK/renormalization_factor_Tmatrix_so.f90 b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
index e9eb7dd..e06ec78 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix_so.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
@@ -40,7 +40,7 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
     do i=nC+1,nO
       cd = 0
       do c=nO+1,nBas-nR
-        do d=nO+1,c-1
+        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
@@ -55,7 +55,7 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
     do a=nO+1,nBas-nR
       kl = 0
       do k=nC+1,nO
-        do l=nC+1,k-1
+        do l=k+1,nO
           kl = kl + 1
           eps = e(p) + e(a) - Omega2(kl)
           Z(p) = Z(p) - rho2(p,a,kl)**2/eps**2
diff --git a/src/QuAcK/self_energy_Tmatrix_diag_so.f90 b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
index c35dedd..3f52f24 100644
--- a/src/QuAcK/self_energy_Tmatrix_diag_so.f90
+++ b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
@@ -44,7 +44,7 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
     do i=nC+1,nO
       cd = 0
       do c=nO+1,nBas-nR
-        do d=nO+1,c-1
+        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
@@ -59,7 +59,7 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
     do a=nO+1,nBas-nR
       kl = 0
       do k=nC+1,nO
-        do l=nC+1,k-1
+        do l=k+1,nO
           kl = kl + 1
           eps = e(p) + e(a) - Omega2(kl)
           SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90
index 5543682..b958b64 100644
--- a/src/QuAcK/soG0T0.f90
+++ b/src/QuAcK/soG0T0.f90
@@ -47,31 +47,30 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
   write(*,*)'************************************************'
   write(*,*)
 
-! Spatial to spin orbitals
+! Define occupied and virtual spaces
 
   nBas2 = 2*nBas
+  nO2   = 2*nO
+  nV2   = 2*nV
+  nC2   = 2*nC 
+  nR2   = 2*nR
+
+! Spatial to spin orbitals
 
   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
+  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), & 
+  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))
 
@@ -79,7 +78,7 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 ! Spinorbital basis
 !----------------------------------------------
 
- ispin = 2
+ ispin = 3
 
 ! Compute linear response