diff --git a/examples/molecule.C2 b/examples/molecule.C2
index 2d592cf..e00bed3 100644
--- a/examples/molecule.C2
+++ b/examples/molecule.C2
@@ -1,5 +1,5 @@
 # nAt nEla nElb nCore nRyd
- 2   7     7   0     0
+ 2   6     6   0     0
 # Znuc   x            y           z
   C     0.           0.         0.
   C     0.           0.         2.0
diff --git a/input/basis b/input/basis
index b2b2293..7a2f348 100644
--- a/input/basis
+++ b/input/basis
@@ -1,30 +1,32 @@
-1   6
-S   8
-  1   1469.0000000              0.0007660        
-  2    220.5000000              0.0058920        
-  3     50.2600000              0.0296710        
-  4     14.2400000              0.1091800        
-  5      4.5810000              0.2827890        
-  6      1.5800000              0.4531230        
-  7      0.5640000              0.2747740        
-  8      0.0734500              0.0097510        
-S   8
-  1   1469.0000000             -0.0001200        
-  2    220.5000000             -0.0009230        
-  3     50.2600000             -0.0046890        
-  4     14.2400000             -0.0176820        
-  5      4.5810000             -0.0489020        
-  6      1.5800000             -0.0960090        
-  7      0.5640000             -0.1363800        
-  8      0.0734500              0.5751020        
+1  6
+S   9
+1         6.665000E+03           6.920000E-04
+2         1.000000E+03           5.329000E-03
+3         2.280000E+02           2.707700E-02
+4         6.471000E+01           1.017180E-01
+5         2.106000E+01           2.747400E-01
+6         7.495000E+00           4.485640E-01
+7         2.797000E+00           2.850740E-01
+8         5.215000E-01           1.520400E-02
+9         1.596000E-01          -3.191000E-03
+S   9
+1         6.665000E+03          -1.460000E-04
+2         1.000000E+03          -1.154000E-03
+3         2.280000E+02          -5.725000E-03
+4         6.471000E+01          -2.331200E-02
+5         2.106000E+01          -6.395500E-02
+6         7.495000E+00          -1.499810E-01
+7         2.797000E+00          -1.272620E-01
+8         5.215000E-01           5.445290E-01
+9         1.596000E-01           5.804960E-01
 S   1
-  1      0.0280500              1.0000000        
-P   3
-  1      1.5340000              0.0227840        
-  2      0.2749000              0.1391070        
-  3      0.0736200              0.5003750        
+1         1.596000E-01           1.000000E+00
+P   4
+1         9.439000E+00           3.810900E-02
+2         2.002000E+00           2.094800E-01
+3         5.456000E-01           5.085570E-01
+4         1.517000E-01           4.688420E-01
 P   1
-  1      0.0240300              1.0000000        
+1         1.517000E-01           1.000000E+00
 D   1
-  1      0.1239000              1.0000000
-
+1         5.500000E-01           1.0000000
diff --git a/input/methods b/input/methods
index 63c03d4..eb2074b 100644
--- a/input/methods
+++ b/input/methods
@@ -1,7 +1,7 @@
 # RHF UHF MOM
-  T   F   F
+  F   T   F
 # MP2 MP3 MP2-F12
-  F   F   F
+  T   F   F
 # CCD CCSD CCSD(T) 
   F   F    F
 # drCCD rCCD lCCD pCCD
@@ -13,7 +13,7 @@
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    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 058d6dd..7f71f46 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,4 @@
 # nAt nEla nElb nCore nRyd
-  1   2     1   0     0
+ 1   5     1   0     0
 # Znuc   x            y           z
-  Li     0.0          0.0         0.0
+  C     0.           0.         0.
diff --git a/input/molecule.xyz b/input/molecule.xyz
index c9a5a65..8162169 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,3 +1,3 @@
   1
 
- Li     0.0000000000    0.0000000000    0.0000000000
+ C      0.0000000000    0.0000000000    0.0000000000
diff --git a/input/options b/input/options
index 348d4b6..daaa7f4 100644
--- a/input/options
+++ b/input/options
@@ -7,9 +7,9 @@
 # spin: singlet triplet TDA
         T       T       F
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
-      256    0.00001 T    5      T   0.00367493 3      
+      256    0.00001 T    5      T   0.0  3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta        COHSEX SOSEX TDA_W G0W GW0 
-         256    0.00001   T    5     T  0.00367493  F      F     F     F   F  
+         256    0.00001   T    5     T  0.0         F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
diff --git a/src/QuAcK/AOtoMO_integral_transform.f90 b/src/QuAcK/AOtoMO_integral_transform.f90
index 17d4119..23e4198 100644
--- a/src/QuAcK/AOtoMO_integral_transform.f90
+++ b/src/QuAcK/AOtoMO_integral_transform.f90
@@ -1,14 +1,17 @@
-subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
+subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
 
-! AO to MO transformation of two-electron integrals
-! Semi-direct O(N^5) algorithm
+! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
+! bra and ket are the spin of (bra|ket)
 
   implicit none
+  include 'parameters.h'
 
 ! Input variables
 
+  integer,intent(in)            :: bra
+  integer,intent(in)            :: ket
   integer,intent(in)            :: nBas
-  double precision,intent(in)   :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas)
+  double precision,intent(in)   :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin)
 
 ! Local variables
 
@@ -20,8 +23,11 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
   double precision,intent(out)  :: ERI_MO_basis(nBas,nBas,nBas,nBas)
 
 ! Memory allocation
+
   allocate(scr(nBas,nBas,nBas,nBas))
 
+! Four-index transform via semi-direct O(N^5) algorithm
+
   scr(:,:,:,:) = 0d0
 
   do l=1,nBas
@@ -29,7 +35,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
       do la=1,nBas
         do nu=1,nBas
           do mu=1,nBas
-            scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l)
+            scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket)
           enddo
         enddo
       enddo
@@ -43,7 +49,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
       do nu=1,nBas
         do i=1,nBas
           do mu=1,nBas
-            ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i)*scr(mu,nu,la,l)
+            ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra)*scr(mu,nu,la,l)
           enddo
         enddo
       enddo
@@ -57,7 +63,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
       do la=1,nBas
         do nu=1,nBas
           do i=1,nBas
-            scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k)
+            scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra)
           enddo
         enddo
       enddo
@@ -71,7 +77,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
       do j=1,nBas
         do i=1,nBas
           do nu=1,nBas
-            ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j)*scr(i,nu,k,l)
+            ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket)*scr(i,nu,k,l)
           enddo
 !         print*,i,k,j,l,ERI_MO_basis(i,j,k,l)
         enddo
diff --git a/src/QuAcK/MP2.f90 b/src/QuAcK/MP2.f90
index 3f48a11..2c6d75c 100644
--- a/src/QuAcK/MP2.f90
+++ b/src/QuAcK/MP2.f90
@@ -1,6 +1,6 @@
 subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
 
-! Perform third-order Moller-Plesset calculation
+! Perform second-order Moller-Plesset calculation
 
   implicit none
 
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 1b3621d..263f819 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -4,6 +4,7 @@ program QuAcK
   include 'parameters.h'
 
   logical                       :: doSph
+  logical                       :: unrestricted
   logical                       :: doRHF,doUHF,doMOM 
   logical                       :: doMP2,doMP3,doMP2F12
   logical                       :: doCCD,doCCSD,doCCSDT
@@ -27,8 +28,8 @@ program QuAcK
   double precision,allocatable  :: ZNuc(:),rNuc(:,:)
   double precision,allocatable  :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
 
-  double precision,allocatable  :: eG0W0(:)
-  double precision,allocatable  :: eG0T0(:)
+  double precision,allocatable  :: eG0W0(:,:)
+  double precision,allocatable  :: eG0T0(:,:)
 
   logical                       :: doACFDT
   logical                       :: exchange_kernel
@@ -50,6 +51,11 @@ program QuAcK
   double precision,allocatable  :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:)
   double precision,allocatable  :: ERI_AO(:,:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
+  integer                       :: bra
+  integer                       :: ket
+  double precision,allocatable  :: ERI_MO_aa(:,:,:,:)
+  double precision,allocatable  :: ERI_MO_ab(:,:,:,:)
+  double precision,allocatable  :: ERI_MO_bb(:,:,:,:)
   double precision,allocatable  :: ERI_ERF_AO(:,:,:,:)
   double precision,allocatable  :: ERI_ERF_MO(:,:,:,:)
   double precision,allocatable  :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
@@ -211,9 +217,8 @@ program QuAcK
 
 ! Memory allocation for one- and two-electron integrals
 
-  allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),eG0T0(nBas),PHF(nBas,nBas,nspin),          &
-           S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), &
-           ERI_AO(nBas,nBas,nBas,nBas),ERI_MO(nBas,nBas,nBas,nBas))
+  allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), &
+           S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas))
 
 ! Read integrals
 
@@ -263,6 +268,9 @@ program QuAcK
 
   if(doUHF) then
 
+    ! Switch on the unrestricted flag
+    unrestricted = .true.
+
     call cpu_time(start_HF)
     call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,EUHF,eHF,cHF,PHF)
     call cpu_time(end_HF)
@@ -307,13 +315,50 @@ program QuAcK
 
   if(doSph) then
 
+    allocate(ERI_MO(nBas,nBas,nBas,nBas))
     ERI_MO(:,:,:,:) = ERI_AO(:,:,:,:)
     print*,'!!! MO = AO !!!'
     deallocate(ERI_AO)
 
   else
 
-    call AOtoMO_integral_transform(nBas,cHF,ERI_AO,ERI_MO)
+    if(unrestricted) then
+
+      ! Memory allocation
+     
+      allocate(ERI_MO_aa(nBas,nBas,nBas,nBas),ERI_MO_ab(nBas,nBas,nBas,nBas),ERI_MO_bb(nBas,nBas,nBas,nBas))
+     
+      ! 4-index transform for (aa|aa) block
+     
+      bra = 1
+      ket = 1
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aa)
+      
+      ! 4-index transform for (bb|bb) block
+     
+      bra = 1
+      ket = 2
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_ab)
+     
+      ! 4-index transform for (aa|bb) block
+     
+      bra = 2
+      ket = 2
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bb)
+     
+    else
+
+      ! Memory allocation
+     
+      allocate(ERI_MO(nBas,nBas,nBas,nBas))
+     
+      ! 4-index transform 
+     
+      bra = 1
+      ket = 1
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO)
+
+    end if
 
   end if
 
@@ -330,7 +375,17 @@ program QuAcK
   if(doMP2) then
 
     call cpu_time(start_MP2)
-    call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
+
+    if(unrestricted) then
+
+      call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,ENuc,EUHF,eHF,EcMP2)
+
+    else
+
+      call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
+
+    end if
+
     call cpu_time(end_MP2)
 
     t_MP2 = end_MP2 - start_MP2
@@ -670,14 +725,24 @@ program QuAcK
 ! Perform G0W0 calculatiom
 !------------------------------------------------------------------------
 
-  eG0W0(:) = eHF(:,1)
+  eG0W0(:,:) = eHF(:,:)
 
   if(doG0W0) then
     
     call cpu_time(start_G0W0)
-    call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 
-              dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & 
-              nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
+    if(unrestricted) then 
+
+      call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,       & 
+                 dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & 
+                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
+    else
+
+      call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 
+                dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & 
+                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
+
+    end if
+
     call cpu_time(end_G0W0)
   
     t_G0W0 = end_G0W0 - start_G0W0
@@ -726,7 +791,7 @@ program QuAcK
 ! Perform G0T0 calculatiom
 !------------------------------------------------------------------------
 
-  eG0T0(:) = eHF(:,1)
+  eG0T0(:,:) = eHF(:,:)
 
   if(doG0T0) then
     
@@ -868,7 +933,7 @@ program QuAcK
     call cpu_time(start_G0W0)
     call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 
               dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & 
-              nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
+              nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
     call cpu_time(end_G0W0)
   
     t_G0W0 = end_G0W0 - start_G0W0
@@ -889,7 +954,7 @@ program QuAcK
     write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
     write(*,*)
 
-    call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
+!   call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
 
   end if
 
diff --git a/src/QuAcK/UMP2.f90 b/src/QuAcK/UMP2.f90
new file mode 100644
index 0000000..674eae9
--- /dev/null
+++ b/src/QuAcK/UMP2.f90
@@ -0,0 +1,158 @@
+subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
+
+! Perform unrestricted second-order Moller-Plesset calculation
+
+  implicit none
+  include 'parameters.h'
+
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: e(nBas,nspin)
+
+! Local variables
+
+  integer                       :: bra,ket
+  integer                       :: i,j,a,b
+  double precision              :: eps
+  double precision              :: Edaa,Exaa,Ecaa
+  double precision              :: Edab,Exab,Ecab
+  double precision              :: Edbb,Exbb,Ecbb
+  double precision              :: Ed,Ex
+
+! Output variables
+
+  double precision,intent(out)  :: Ec
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'********************************************************'
+  write(*,*)'| Unrestricted second-order Moller-Plesset calculation |'
+  write(*,*)'********************************************************'
+  write(*,*)
+
+!---------------------!
+! Compute UMP2 energy |
+!---------------------!
+
+! aaaa block
+
+  bra = 1
+  ket = 1
+
+  Edaa = 0d0
+  Exaa = 0d0
+
+  do i=nC(bra)+1,nO(bra)
+    do a=nO(bra)+1,nBas-nR(bra)
+
+      do j=nC(ket)+1,nO(ket)
+        do b=nO(ket)+1,nBas-nR(ket)
+
+          eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket) 
+         
+          Edaa = Edaa + 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,a,b)/eps
+          Exaa = Exaa - 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,b,a)/eps
+
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  Ecaa = Edaa + Exaa
+
+! aabb block
+
+  bra = 1
+  ket = 2
+
+  Edab = 0d0
+  Exab = 0d0
+
+  do i=nC(bra)+1,nO(bra)
+    do a=nO(bra)+1,nBas-nR(bra)
+
+      do j=nC(ket)+1,nO(ket)
+        do b=nO(ket)+1,nBas-nR(ket)
+
+          eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket) 
+         
+          Edab = Edab + ERI_ab(i,j,a,b)*ERI_ab(i,j,a,b)/eps
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  Ecab = Edab + Exab
+
+! bbbb block
+
+  bra = 2
+  ket = 2
+
+  Edbb = 0d0
+  Exbb = 0d0
+
+  do i=nC(bra)+1,nO(bra)
+    do a=nO(bra)+1,nBas-nR(bra)
+
+      do j=nC(ket)+1,nO(ket)
+        do b=nO(ket)+1,nBas-nR(ket)
+
+          eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket) 
+         
+          Edbb = Edbb + 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,a,b)/eps
+          Exbb = Exbb - 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,b,a)/eps
+
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  Ecbb = Edbb + Exbb
+
+! Final flush
+
+  Ed = Edaa + Edab + Edbb
+  Ex = Exaa + Exab + Exbb
+  Ec = Ed + Ex
+
+  write(*,*)
+  write(*,'(A32)')           '--------------------------'
+  write(*,'(A32)')           ' MP2 calculation          '
+  write(*,'(A32)')           '--------------------------'
+  write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ec
+  write(*,'(A32,1X,F16.10)') '   alpha-alpha          = ',Ecaa
+  write(*,'(A32,1X,F16.10)') '   alpha-beta           = ',Ecab
+  write(*,'(A32,1X,F16.10)') '    beta-beta           = ',Ecbb
+  write(*,*)
+  write(*,'(A32,1X,F16.10)') ' Direct part            = ',Ed
+  write(*,'(A32,1X,F16.10)') '   alpha-alpha          = ',Edaa
+  write(*,'(A32,1X,F16.10)') '   alpha-beta           = ',Edab
+  write(*,'(A32,1X,F16.10)') '    beta-beta           = ',Edbb
+  write(*,*)
+  write(*,'(A32,1X,F16.10)') ' Exchange part          = ',Ex
+  write(*,'(A32,1X,F16.10)') '   alpha-alpha          = ',Exaa
+  write(*,'(A32,1X,F16.10)') '   alpha-beta           = ',Exab
+  write(*,'(A32,1X,F16.10)') '    beta-beta           = ',Exbb
+  write(*,'(A32)')           '--------------------------'
+  write(*,'(A32,1X,F16.10)') ' MP2 electronic  energy = ',       EHF + Ec
+  write(*,'(A32,1X,F16.10)') ' MP2 total       energy = ',ENuc + EHF + Ec
+  write(*,'(A32)')           '--------------------------'
+  write(*,*)
+
+end subroutine UMP2
diff --git a/src/QuAcK/linear_response_A_matrix.f90 b/src/QuAcK/linear_response_A_matrix.f90
index f2c20e0..efc6e81 100644
--- a/src/QuAcK/linear_response_A_matrix.f90
+++ b/src/QuAcK/linear_response_A_matrix.f90
@@ -15,7 +15,7 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
   
 ! Local variables
 
-  double precision              :: delta_spin,delta_dRPA
+  double precision              :: delta_dRPA
   double precision,external     :: Kronecker_delta
 
   integer                       :: i,j,a,b,ia,jb
@@ -24,35 +24,78 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
 
   double precision,intent(out)  :: A_lr(nS,nS)
 
-! Singlet or triplet manifold?
-
-  delta_spin = 0d0
-  if(ispin == 1) delta_spin = +1d0
-  if(ispin == 2) delta_spin = -1d0
-
 ! Direct RPA
 
   delta_dRPA = 0d0
   if(dRPA) delta_dRPA = 1d0
 
-! Build A matrix
+! Build A matrix for single manifold
 
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
+  if(ispin == 1) then 
 
-          A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                      + (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
-                      - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        ia = ia + 1
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+ 
+            A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                        + 2d0*lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
 
-        enddo
-      enddo
-    enddo
-  enddo
+          end  do
+        end  do
+      end  do
+    end  do
+
+  end if
+
+! Build A matrix for triplet manifold
+
+  if(ispin == 2) then 
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        ia = ia + 1
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+ 
+            A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                        - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+
+          end  do
+        end  do
+      end  do
+    end  do
+
+  end if
+
+! Build A matrix for spin orbitals
+
+  if(ispin == 3) then 
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        ia = ia + 1
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+ 
+            A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                        + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+
+          end  do
+        end  do
+      end  do
+    end  do
+
+  end if
 
 end subroutine linear_response_A_matrix
diff --git a/src/QuAcK/linear_response_B_matrix.f90 b/src/QuAcK/linear_response_B_matrix.f90
index 5226059..5f15159 100644
--- a/src/QuAcK/linear_response_B_matrix.f90
+++ b/src/QuAcK/linear_response_B_matrix.f90
@@ -14,7 +14,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
   
 ! Local variables
 
-  double precision              :: delta_spin,delta_dRPA
+  double precision              :: delta_dRPA
 
   integer                       :: i,j,a,b,ia,jb
 
@@ -22,34 +22,75 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
 
   double precision,intent(out)  :: B_lr(nS,nS)
 
-! Singlet or triplet manifold?
-
-  delta_spin = 0d0
-  if(ispin == 1) delta_spin = +1d0
-  if(ispin == 2) delta_spin = -1d0
-
 ! Direct RPA
 
   delta_dRPA = 0d0
   if(dRPA) delta_dRPA = 1d0
 
-! Build B matrix
+! Build B matrix for singlet manifold
 
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
+  if(ispin == 1) then
 
-          B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
-                      - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        ia = ia + 1
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+ 
+            B_lr(ia,jb) = 2d0*lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+ 
+          end do
+        end do
+      end do
+    end do
 
-        enddo
-      enddo
-    enddo
-  enddo
+  end if
+
+! Build B matrix for triplet manifold
+
+  if(ispin == 2) then
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        ia = ia + 1
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+ 
+            B_lr(ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
+
+! Build B matrix for spin orbitals 
+
+  if(ispin == 3) then
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        ia = ia + 1
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+ 
+            B_lr(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
 
 end subroutine linear_response_B_matrix