diff --git a/input/options b/input/options
index b5853ef..ceeb2e0 100644
--- a/input/options
+++ b/input/options
@@ -15,6 +15,6 @@
 # ACFDT: AC Kx  XBS
          F  F   F
 # BSE: BSE dBSE dTDA evDyn
-        F    F    T    F
+        T    F    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90
index d94379d..1a0e2b0 100644
--- a/src/GT/Bethe_Salpeter_Tmatrix.f90
+++ b/src/GT/Bethe_Salpeter_Tmatrix.f90
@@ -88,10 +88,10 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
                call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TAs)
   if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TBs)
 
-  print*,'ab block of TAxs'
-  call matout(nS,nS,TAs)
-  print*,'ab block of TB'
-  call matout(nS,nS,TBs)
+! print*,'ab block of TA'
+! call matout(nS,nS,TAs)
+! print*,'ab block of TB'
+! call matout(nS,nS,TBs)
 
 !----------------------------------------------
 ! Compute T-matrix for alpha-alpha block
@@ -108,16 +108,10 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
                call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TAt)
   if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TBt)
 
-  print*,'aa block of TA'
-  call matout(nS,nS,TAt)
-  print*,'aa block of TB'
-  call matout(nS,nS,TBt)
-
-  TAs(:,:) = TAs(:,:) - TAt(:,:)
-  TBs(:,:) = TBs(:,:) - TBt(:,:)
-
-  TAt(:,:) = - TAs(:,:)
-  TBt(:,:) = - TBs(:,:)
+! print*,'aa block of TA'
+! call matout(nS,nS,TAt)
+! print*,'aa block of TB'
+! call matout(nS,nS,TBt)
 
 !-------------------
 ! Singlet manifold
@@ -130,7 +124,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
     ! Compute BSE singlet excitation energies
 
-    call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAs,TBs, &
+    call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAs+TAt,TBs+TBt, &
                                  EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
 
     call print_excitation('BSE@GT      ',ispin,nS,OmBSE(:,ispin))
@@ -171,7 +165,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
     ! Compute BSE triplet excitation energies
 
-    call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt,TBt, &
+    call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt-TAs,TBt-TBs, &
                                  EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
     call print_excitation('BSE@GT      ',ispin,nS,OmBSE(:,ispin))
     call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & 
diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90
index a472c1c..0b84f93 100644
--- a/src/GT/G0T0.f90
+++ b/src/GT/G0T0.f90
@@ -80,10 +80,10 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 
 ! Dimensions of the pp-RPA linear reponse matrices
 
-  nOOs = nO*(nO + 1)/2
-  nVVs = nV*(nV + 1)/2
-! nOOs = nO*nO
-! nVVs = nV*nV
+! nOOs = nO*(nO + 1)/2
+! nVVs = nV*(nV + 1)/2
+  nOOs = nO*nO
+  nVVs = nV*nV
 
   nOOt = nO*(nO - 1)/2
   nVVt = nV*(nV - 1)/2
@@ -103,8 +103,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 !----------------------------------------------
 
   ispin  = 1
-  iblock = 1
-! iblock = 3
+! iblock = 1
+  iblock = 3
 
 ! Compute linear response
 
@@ -142,8 +142,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   SigT(:) = 0d0
   Z(:)    = 0d0
 
-  iblock = 1
-! iblock = 3
+! iblock = 1
+  iblock = 3
 
   call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
@@ -189,8 +189,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 ! Compute the ppRPA correlation energy
 
   ispin  = 1
-! iblock = 3
-  iblock = 1
+  iblock = 3
+! iblock = 1
 
   call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eG0T0,ERI_MO,  & 
                           Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
diff --git a/src/LR/linear_response_Tmatrix.f90 b/src/LR/linear_response_Tmatrix.f90
index e2a7156..dbb0596 100644
--- a/src/LR/linear_response_Tmatrix.f90
+++ b/src/LR/linear_response_Tmatrix.f90
@@ -48,8 +48,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
   call matout(nS,nS,A_BSE)
 
   A(:,:) = A(:,:) + A_BSE(:,:)
-! if(ispin == 1) A(:,:) = A(:,:) + A_BSE(:,:)
-! if(ispin == 2) A(:,:) = A(:,:) - A_BSE(:,:)
 
 ! Tamm-Dancoff approximation
 
@@ -71,8 +69,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
     call matout(nS,nS,B_BSE)
 
     B(:,:) = B(:,:) + B_BSE(:,:)
-!   if(ispin == 1) B(:,:) = B(:,:) + B_BSE(:,:)
-!   if(ispin == 2) B(:,:) = B(:,:) - B_BSE(:,:)
 
     ! Build A + B and A - B matrices