diff --git a/input/methods b/input/methods
index 0b7bf09..df12e6b 100644
--- a/input/methods
+++ b/input/methods
@@ -1,7 +1,7 @@
 # RHF UHF KS MOM 
   T   F   F  F   
 # MP2* MP3 MP2-F12
-  F   F   F
+  T   F   F
 # CCD pCCD DCD CCSD CCSD(T) 
   F   F    F   F    F
 # drCCD rCCD crCCD lCCD
@@ -15,7 +15,7 @@
 # G0W0* evGW* qsGW* ufG0W0 ufGW
   F     F     F     F      F
 # G0T0 evGT qsGT 
-  F    F    F
+  T    F    F
 # MCMP2
   F
 # * unrestricted version available
diff --git a/input/options b/input/options
index ceeb2e0..b5853ef 100644
--- a/input/options
+++ b/input/options
@@ -15,6 +15,6 @@
 # ACFDT: AC Kx  XBS
          F  F   F
 # BSE: BSE dBSE dTDA evDyn
-        T    F    T    F
+        F    F    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/mol/h2.xyz b/mol/h2.xyz
index 8832854..d955cc4 100644
--- a/mol/h2.xyz
+++ b/mol/h2.xyz
@@ -1,4 +1,4 @@
 2
  
 H 0. 0. 0.
-H 0. 0. 1.14
+H 0. 0. 0.7
diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90
index b314eb6..d94379d 100644
--- a/src/GT/Bethe_Salpeter_Tmatrix.f90
+++ b/src/GT/Bethe_Salpeter_Tmatrix.f90
@@ -58,7 +58,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
   integer                       :: iblock
 
   double precision              :: EcRPA(nspin)
-  double precision,allocatable  :: TA(:,:),TB(:,:)
+  double precision,allocatable  :: TAs(:,:),TBs(:,:)
+  double precision,allocatable  :: TAt(:,:),TBt(:,:)
   double precision,allocatable  :: OmBSE(:,:)
   double precision,allocatable  :: XpY_BSE(:,:,:)
   double precision,allocatable  :: XmY_BSE(:,:,:)
@@ -69,12 +70,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
 ! Memory allocation
 
-  allocate(TA(nS,nS),TB(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin))
-
-! Initialize T matrix
-
-  TA(:,:) = 0d0
-  TB(:,:) = 0d0
+  allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), & 
+           OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin))
 
 !----------------------------------------------
 ! Compute T-matrix for alpha-beta block
@@ -88,13 +85,13 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
 ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-               call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
-  if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
+               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*,'aa block of TA'
-  call matout(nS,nS,TA)
-  print*,'aa block of TB'
-  call matout(nS,nS,TB)
+  print*,'ab block of TAxs'
+  call matout(nS,nS,TAs)
+  print*,'ab block of TB'
+  call matout(nS,nS,TBs)
 
 !----------------------------------------------
 ! Compute T-matrix for alpha-alpha block
@@ -108,13 +105,19 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
 ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-               call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
-  if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
+               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+ab block of TA'
-  call matout(nS,nS,TA)
-  print*,'aa+ab block of TB'
-  call matout(nS,nS,TB)
+  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(:,:)
 
 !-------------------
 ! Singlet manifold
@@ -127,7 +130,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,TA,TB, &
+    call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAs,TBs, &
                                  EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
 
     call print_excitation('BSE@GT      ',ispin,nS,OmBSE(:,ispin))
@@ -168,7 +171,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,TA,TB, &
+    call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt,TBt, &
                                  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 16225b3..03db7ea 100644
--- a/src/GT/G0T0.f90
+++ b/src/GT/G0T0.f90
@@ -80,6 +80,8 @@ 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
 
@@ -101,6 +103,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 !----------------------------------------------
 
   ispin  = 1
+! iblock = 1
   iblock = 3
 
 ! Compute linear response
@@ -139,7 +142,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   SigT(:) = 0d0
   Z(:)    = 0d0
 
-  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)
 
@@ -147,7 +151,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 
   call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z)
 
-  iblock =  4
+  iblock = 4
 
   call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
@@ -186,6 +190,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 
   ispin  = 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))
   ispin  = 2
diff --git a/src/GT/excitation_density_Tmatrix.f90 b/src/GT/excitation_density_Tmatrix.f90
index 01ae36c..33c0ff8 100644
--- a/src/GT/excitation_density_Tmatrix.f90
+++ b/src/GT/excitation_density_Tmatrix.f90
@@ -53,7 +53,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
           do c=nO+1,nBas-nR
             do d=c,nBas-nR
               cd = cd + 1
-              rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab) 
+              rho1(p,q,ab) = rho1(p,q,ab) & 
+                           + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d)))
             end do
           end do
  
@@ -61,7 +62,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
           do k=nC+1,nO
             do l=k,nO
               kl = kl + 1
-              rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab) 
+              rho1(p,q,ab) = rho1(p,q,ab) & 
+                           + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l)))
             end do
           end do
  
@@ -73,7 +75,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
           do c=nO+1,nBas-nR
             do d=c,nBas-nR
               cd = cd + 1
-              rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij)
+              rho2(p,q,ij) = rho2(p,q,ij) &
+                           + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d)))
             end do
           end do
  
@@ -81,7 +84,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
           do k=nC+1,nO
             do l=k,nO
               kl = kl + 1
-              rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij) 
+              rho2(p,q,ij) = rho2(p,q,ij) &
+                           + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l)))
             end do
           end do
  
diff --git a/src/GT/self_energy_Tmatrix.f90 b/src/GT/self_energy_Tmatrix.f90
index b5d82a4..13c211f 100644
--- a/src/GT/self_energy_Tmatrix.f90
+++ b/src/GT/self_energy_Tmatrix.f90
@@ -69,7 +69,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
     do j=nC+1,nO
       do cd=1,nVV
         eps = e(i) + e(j) - Omega1(cd)
-        EcGM = EcGM - rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
+        EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
       enddo
     enddo
   enddo
@@ -78,7 +78,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
     do b=nO+1,nBas-nR
       do kl=1,nOO
         eps = e(a) + e(b) - Omega2(kl)
-        EcGM = EcGM + rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+        EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
       enddo
     enddo
   enddo
diff --git a/src/GT/self_energy_Tmatrix_diag.f90 b/src/GT/self_energy_Tmatrix_diag.f90
index 39a3bfc..c52eb2d 100644
--- a/src/GT/self_energy_Tmatrix_diag.f90
+++ b/src/GT/self_energy_Tmatrix_diag.f90
@@ -65,7 +65,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
     do j=nC+1,nO
       do cd=1,nVV
         eps = e(i) + e(j) - Omega1(cd)
-        EcGM = EcGM - rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
+        EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
       enddo
     enddo
   enddo
@@ -74,7 +74,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
     do b=nO+1,nBas-nR
       do kl=1,nOO
         eps = e(a) + e(b) - Omega2(kl)
-        EcGM = EcGM + rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+        EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
       enddo
     enddo
   enddo
diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90
index 9e3505e..2e1586a 100644
--- a/src/GT/static_Tmatrix_A.f90
+++ b/src/GT/static_Tmatrix_A.f90
@@ -34,6 +34,8 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome
 
   double precision,intent(out)  :: TA(nS,nS)
 
+  TA(:,:) = 0d0
+
   ia = 0
   do i=nC+1,nO
     do a=nO+1,nBas-nR
diff --git a/src/GT/static_Tmatrix_B.f90 b/src/GT/static_Tmatrix_B.f90
index 7cf9bb0..01ef3ab 100644
--- a/src/GT/static_Tmatrix_B.f90
+++ b/src/GT/static_Tmatrix_B.f90
@@ -34,6 +34,8 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome
 
   double precision,intent(out)  :: TB(nS,nS)
 
+  TB(:,:) = 0d0
+
   ia = 0
   do i=nC+1,nO
     do a=nO+1,nBas-nR
diff --git a/src/LR/linear_response_Tmatrix.f90 b/src/LR/linear_response_Tmatrix.f90
index 72d3750..e2a7156 100644
--- a/src/LR/linear_response_Tmatrix.f90
+++ b/src/LR/linear_response_Tmatrix.f90
@@ -42,15 +42,15 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
 
   call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
 
+  print*,'A'
+  call matout(nS,nS,A)
+  print*,'TA'
+  call matout(nS,nS,A_BSE)
+
   A(:,:) = A(:,:) + A_BSE(:,:)
 ! if(ispin == 1) A(:,:) = A(:,:) + A_BSE(:,:)
 ! if(ispin == 2) A(:,:) = A(:,:) - A_BSE(:,:)
 
-! print*,'A'
-! call matout(nS,nS,A)
-! print*,'TA'
-! call matout(nS,nS,A_BSE)
-
 ! Tamm-Dancoff approximation
 
   if(TDA) then
@@ -65,15 +65,15 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
 
     call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
 
+    print*,'B'
+    call matout(nS,nS,B)
+    print*,'TB'
+    call matout(nS,nS,B_BSE)
+
     B(:,:) = B(:,:) + B_BSE(:,:)
 !   if(ispin == 1) B(:,:) = B(:,:) + B_BSE(:,:)
 !   if(ispin == 2) B(:,:) = B(:,:) - B_BSE(:,:)
 
-! print*,'B'
-! call matout(nS,nS,B)
-! print*,'TB'
-! call matout(nS,nS,B_BSE)
-
     ! Build A + B and A - B matrices 
 
     ApB = A + B