diff --git a/input/methods b/input/methods
index 762c3e7..ff3db34 100644
--- a/input/methods
+++ b/input/methods
@@ -15,7 +15,7 @@
 # G0W0* evGW* qsGW*
   F    F    F
 # G0T0 evGT qsGT 
-  T    F    F
+  T    T    T
 # MCMP2
   F
 # * unrestricted version available
diff --git a/input/options b/input/options
index 51eb71f..6de1fec 100644
--- a/input/options
+++ b/input/options
@@ -13,6 +13,6 @@
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    T    T    F
+        F    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 b/src/MBPT/Bethe_Salpeter_Tmatrix.f90
index 94541b3..07a192f 100644
--- a/src/MBPT/Bethe_Salpeter_Tmatrix.f90
+++ b/src/MBPT/Bethe_Salpeter_Tmatrix.f90
@@ -56,8 +56,6 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
   integer                       :: ispin
   integer                       :: iblock
-  integer                       :: dERI
-  integer                       :: xERI
 
   double precision              :: EcRPA(nspin)
   double precision,allocatable  :: TA(:,:),TB(:,:)
@@ -84,13 +82,11 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
   ispin  = 1
   iblock = 3
-  dERI   = +1d0
-  xERI   = +0d0
 
   call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eT,ERI,  &
                           Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
 
-! call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
                call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
   if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
@@ -101,13 +97,11 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
 
   ispin  = 2
   iblock = 4
-  dERI   = +1d0
-  xERI   = -1d0
 
   call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eT,ERI,  &
                           Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
 
-! call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
                call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
   if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90
index 3574a56..06d1b76 100644
--- a/src/MBPT/G0T0.f90
+++ b/src/MBPT/G0T0.f90
@@ -44,9 +44,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   integer                       :: iblock
   integer                       :: nOOs,nOOt
   integer                       :: nVVs,nVVt
-  double precision              :: dERI
-  double precision              :: xERI
-  double precision              :: alpha
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
   double precision              :: EcAC(nspin)
@@ -140,20 +137,16 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   Z(:)    = 0d0
 
   iblock =  3
-  dERI   = +1d0
-  xERI   = +0d0
 
-  call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+  call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
   call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,SigT)
 
   call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z)
 
   iblock =  4
-  dERI   = +1d0
-  xERI   = -1d0
 
-  call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+  call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
   call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,SigT)
 
diff --git a/src/MBPT/evGT.f90 b/src/MBPT/evGT.f90
index 64e3b41..eaa8c35 100644
--- a/src/MBPT/evGT.f90
+++ b/src/MBPT/evGT.f90
@@ -54,9 +54,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
   integer                       :: iblock
   integer                       :: nOOs,nOOt
   integer                       :: nVVs,nVVt
-  double precision              :: dERI
-  double precision              :: xERI
-  double precision              :: alpha
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
   double precision              :: EcAC(nspin)
@@ -76,11 +73,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
   double precision,allocatable  :: SigT(:)
   double precision,allocatable  :: Z(:)
 
-  double precision,allocatable  :: Omega(:,:)
-  double precision,allocatable  :: XpY(:,:,:)
-  double precision,allocatable  :: XmY(:,:,:)
-  double precision,allocatable  :: rho(:,:,:,:)
-
 ! Output variables
 
 ! Hello world
@@ -140,14 +132,9 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
  
     ! Compute linear response
  
-    call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:),ERI_MO(:,:,:,:),  & 
-                            Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin))
+    call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO,  & 
+                            Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
  
-!   EcRPA(ispin) = 1d0*EcRPA(ispin)
- 
-!   call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:))
-!   call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:))
-
   !----------------------------------------------
   ! alpha-alpha block
   !----------------------------------------------
@@ -157,14 +144,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
 
   ! Compute linear response
 
-    call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:),ERI_MO(:,:,:,:),  & 
-                          Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin))
-
-!   EcRPA(ispin) = 2d0*EcRPA(ispin)
-!   EcRPA(ispin) = 3d0*EcRPA(ispin)
-
-!   call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:))
-!   call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:))
+    call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO,  & 
+                          Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
 
   !----------------------------------------------
   ! Compute T-matrix version of the self-energy 
@@ -174,32 +155,26 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
     Z(:)    = 0d0
  
     iblock =  3
-    dERI   = +1d0
-    xERI   = +0d0
-    alpha  = +1d0
  
-    call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO(:,:,:,:), &
-                                    X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:))
+    call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, &
+                                    X1s,Y1s,rho1s,X2s,Y2s,rho2s)
  
-    call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:), & 
-                                  Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),SigT(:))
+    call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & 
+                                  Omega1s,rho1s,Omega2s,rho2s,SigT)
  
-    call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:), & 
-                                        Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),Z(:))
+    call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & 
+                                        Omega1s,rho1s,Omega2s,rho2s,Z)
  
     iblock =  4
-    dERI   = +1d0
-    xERI   = -1d0
-    alpha  = +1d0
  
-    call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO(:,:,:,:), &
-                                    X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:))
+    call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, &
+                                    X1t,Y1t,rho1t,X2t,Y2t,rho2t)
  
-    call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), & 
-                                  Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:))
+    call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & 
+                                  Omega1t,rho1t,Omega2t,rho2t,SigT)
  
-    call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), & 
-                                        Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),Z(:))
+    call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & 
+                                        Omega1t,rho1t,Omega2t,rho2t,Z)
  
     Z(:) = 1d0/(1d0 - Z(:))
 
@@ -219,7 +194,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
   ! Dump results
   !----------------------------------------------
 
-    call print_evGT(nBas,nO,nSCF,Conv,eHF(:),SigT(:),Z(:),eGT(:))
+    call print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
 
     ! DIIS extrapolation
 
@@ -247,12 +222,12 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
 
   ispin  = 1
   iblock = 3
-  call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:),ERI_MO(:,:,:,:),  &
-                          Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin))
+  call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO,  &
+                          Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
   ispin  = 2
   iblock = 4
-  call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:),ERI_MO(:,:,:,:),  &
-                          Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin))
+  call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO,  &
+                          Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
   EcRPA(1) = EcRPA(1) - EcRPA(2)
   EcRPA(2) = 3d0*EcRPA(2)
 
diff --git a/src/MBPT/excitation_density_Tmatrix.f90 b/src/MBPT/excitation_density_Tmatrix.f90
index 9079e10..01ae36c 100644
--- a/src/MBPT/excitation_density_Tmatrix.f90
+++ b/src/MBPT/excitation_density_Tmatrix.f90
@@ -1,4 +1,4 @@
-subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
+subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
 
 ! Compute excitation densities for T-matrix self-energy
 
@@ -7,8 +7,6 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
 ! Input variables
 
   integer,intent(in)            :: ispin
-  double precision,intent(in)   :: dERI
-  double precision,intent(in)   :: xERI
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -55,8 +53,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do c=nO+1,nBas-nR
             do d=c,nBas-nR
               cd = cd + 1
-              rho1(p,q,ab) = rho1(p,q,ab) & 
-                           + (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X1(cd,ab) 
+              rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab) 
             end do
           end do
  
@@ -64,8 +61,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do k=nC+1,nO
             do l=k,nO
               kl = kl + 1
-              rho1(p,q,ab) = rho1(p,q,ab) & 
-                           + (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y1(kl,ab) 
+              rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab) 
             end do
           end do
  
@@ -77,8 +73,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do c=nO+1,nBas-nR
             do d=c,nBas-nR
               cd = cd + 1
-              rho2(p,q,ij) = rho2(p,q,ij) & 
-                           + (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X2(cd,ij)
+              rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij)
             end do
           end do
  
@@ -86,8 +81,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do k=nC+1,nO
             do l=k,nO
               kl = kl + 1
-              rho2(p,q,ij) = rho2(p,q,ij) & 
-                           + (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y2(kl,ij) 
+              rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij) 
             end do
           end do
  
@@ -114,7 +108,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
             do d=c+1,nBas-nR
               cd = cd + 1
               rho1(p,q,ab) = rho1(p,q,ab) & 
-                           + (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X1(cd,ab) 
+                           + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) 
             end do
           end do
  
@@ -123,7 +117,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
             do l=k+1,nO
               kl = kl + 1
               rho1(p,q,ab) = rho1(p,q,ab) & 
-                           + (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y1(kl,ab) 
+                           + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) 
             end do
           end do
  
@@ -136,7 +130,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
             do d=c+1,nBas-nR
               cd = cd + 1
               rho2(p,q,ij) = rho2(p,q,ij) & 
-                           + (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X2(cd,ij) 
+                           + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) 
             end do
           end do
  
@@ -145,7 +139,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
             do l=k+1,nO
               kl = kl + 1
               rho2(p,q,ij) = rho2(p,q,ij) & 
-                           + (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y2(kl,ij) 
+                           + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) 
             end do
           end do
  
@@ -171,8 +165,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do c=nO+1,nBas-nR
             do d=nO+1,nBas-nR
               cd = cd + 1
-              rho1(p,q,ab) = rho1(p,q,ab) & 
-                           + (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X1(cd,ab) 
+              rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab) 
             end do
           end do
  
@@ -180,8 +173,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do k=nC+1,nO
             do l=nC+1,nO
               kl = kl + 1
-              rho1(p,q,ab) = rho1(p,q,ab) & 
-                           + (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y1(kl,ab) 
+              rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab) 
             end do
           end do
  
@@ -193,8 +185,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do c=nO+1,nBas-nR
             do d=nO+1,nBas-nR
               cd = cd + 1
-              rho2(p,q,ij) = rho2(p,q,ij) & 
-                           + (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X2(cd,ij) 
+              rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij) 
             end do
           end do
  
@@ -202,8 +193,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
           do k=nC+1,nO
             do l=nC+1,nO
               kl = kl + 1
-              rho2(p,q,ij) = rho2(p,q,ij) & 
-                           + (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y2(kl,ij) 
+              rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij) 
             end do
           end do
  
diff --git a/src/MBPT/print_qsGT.f90 b/src/MBPT/print_qsGT.f90
new file mode 100644
index 0000000..5883b62
--- /dev/null
+++ b/src/MBPT/print_qsGT.f90
@@ -0,0 +1,120 @@
+subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
+
+! Print one-electron energies and other stuff for qsGT
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
+  integer,intent(in)                 :: nSCF
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ET
+  double precision,intent(in)        :: EV
+  double precision,intent(in)        :: EJ
+  double precision,intent(in)        :: Ex
+  double precision,intent(in)        :: EcGM
+  double precision,intent(in)        :: EcRPA(nspin)
+  double precision,intent(in)        :: Conv
+  double precision,intent(in)        :: thresh
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: eGT(nBas)
+  double precision,intent(in)        :: c(nBas)
+  double precision,intent(in)        :: SigC(nBas,nBas)
+  double precision,intent(in)        :: Z(nBas)
+  double precision,intent(in)        :: dipole(ncart)
+
+! Local variables
+
+  integer                            :: x,ixyz,HOMO,LUMO
+  double precision                   :: Gap
+  double precision,external          :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)       :: EqsGT
+
+! HOMO and LUMO
+
+  HOMO = nO
+  LUMO = HOMO + 1
+  Gap = eGT(LUMO)-eGT(HOMO)
+
+! Compute energies
+
+! Dump results
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  if(nSCF < 10) then
+    write(*,'(1X,A21,I1,A1,I1,A12)')'  Self-consistent qsG',nSCF,'W',nSCF,' calculation'
+  else
+    write(*,'(1X,A21,I2,A1,I2,A12)')'  Self-consistent qsG',nSCF,'W',nSCF,' calculation'
+  endif
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','Sig_T (eV)','|','Z','|','e_QP (eV)','|'
+  write(*,*)'-------------------------------------------------------------------------------'
+
+  do x=1,nBas
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+    '|',x,'|',eHF(x)*HaToeV,'|',SigC(x,x)*HaToeV,'|',Z(x),'|',eGT(x)*HaToeV,'|'
+  enddo
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A10,I3)')   'Iteration ',nSCF
+  write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
+  write(*,*)'-------------------------------------------'
+  write(*,'(2X,A30,F15.6,A3)') 'qsGT HOMO      energy:',eGT(HOMO)*HaToeV,' eV'
+  write(*,'(2X,A30,F15.6,A3)') 'qsGT LUMO      energy:',eGT(LUMO)*HaToeV,' eV'
+  write(*,'(2X,A30,F15.6,A3)') 'qsGT HOMO-LUMO gap   :',Gap*HaToeV,' eV'
+  write(*,*)'-------------------------------------------'
+  write(*,'(2X,A30,F15.6,A3)') '    qsGT total       energy:',ENuc + EqsGT,' au'
+  write(*,'(2X,A30,F15.6,A3)') '    qsGT exchange    energy:',Ex,' au'
+  write(*,'(2X,A30,F15.6,A3)') ' GM@qsGT correlation energy:',EcGM,' au'
+  write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGT correlation energy:',sum(EcRPA(:)),' au'
+  write(*,*)'-------------------------------------------'
+  write(*,*)
+
+! Dump results for final iteration
+
+  if(Conv < thresh) then
+
+    write(*,*)
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32)')           ' Summary              '
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Kinetic      energy: ',ET,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Potential    energy: ',EV,' au'
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Hartree      energy: ',EJ,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Exchange     energy: ',Ex,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Correlation  energy: ',EcGM,' au'
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32,1X,F16.10,A3)') ' Electronic   energy: ',EqsGT,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Nuclear   repulsion: ',ENuc,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' qsGT         energy: ',ENuc + EqsGT,' au'
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A35)')           ' Dipole moment (Debye)    '
+    write(*,'(10X,4A10)')      'X','Y','Z','Tot.'
+    write(*,'(10X,4F10.6)')    (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
+    write(*,'(A50)')           '-----------------------------------------'
+    write(*,*)
+ 
+    write(*,'(A50)') '---------------------------------------'
+    write(*,'(A32)') ' qsGT MO coefficients'
+    write(*,'(A50)') '---------------------------------------'
+    call matout(nBas,nBas,c)
+    write(*,*)
+    write(*,'(A50)') '---------------------------------------'
+    write(*,'(A32)') ' qsGT MO energies'
+    write(*,'(A50)') '---------------------------------------'
+    call matout(nBas,1,eGT)
+    write(*,*)
+
+  endif
+
+end subroutine print_qsGT
diff --git a/src/MBPT/qsGT.f90 b/src/MBPT/qsGT.f90
new file mode 100644
index 0000000..5f977d6
--- /dev/null
+++ b/src/MBPT/qsGT.f90
@@ -0,0 +1,399 @@
+subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,          & 
+                dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
+                S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+
+! Perform a quasiparticle self-consistent GT calculation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: maxSCF
+  integer,intent(in)            :: max_diis
+  double precision,intent(in)   :: thresh
+  logical,intent(in)            :: doACFDT
+  logical,intent(in)            :: exchange_kernel
+  logical,intent(in)            :: doXBS
+  logical,intent(in)            :: BSE
+  logical,intent(in)            :: TDA_T
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: dBSE
+  logical,intent(in)            :: dTDA
+  logical,intent(in)            :: evDyn
+  logical,intent(in)            :: singlet
+  logical,intent(in)            :: triplet
+  double precision,intent(in)   :: eta
+
+  integer,intent(in)            :: nNuc
+  double precision,intent(in)   :: ZNuc(nNuc)
+  double precision,intent(in)   :: rNuc(nNuc,ncart)
+  double precision,intent(in)   :: ENuc
+
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: cHF(nBas,nBas)
+  double precision,intent(in)   :: PHF(nBas,nBas)
+  double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: Hc(nBas,nBas)
+  double precision,intent(in)   :: X(nBas,nBas)
+  double precision,intent(in)   :: ERI_AO(nBas,nBas,nBas,nBas)
+  double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int_AO(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int_MO(nBas,nBas,ncart)
+
+! Local variables
+
+  integer                       :: nSCF
+  integer                       :: nBasSq
+  integer                       :: ispin
+  integer                       :: iblock
+  integer                       :: n_diis
+  double precision              :: ET
+  double precision              :: EV
+  double precision              :: EJ
+  double precision              :: Ex
+  double precision              :: EqsGT
+  double precision              :: EcRPA(nspin)
+  double precision              :: EcBSE(nspin)
+  double precision              :: EcAC(nspin)
+  double precision              :: EcGM
+  double precision              :: Conv
+  double precision              :: rcond
+  double precision,external     :: trace_matrix
+  double precision              :: dipole(ncart)
+
+  integer                       :: nOOs,nOOt
+  integer                       :: nVVs,nVVt
+
+  logical                       :: print_W = .false.
+  double precision,allocatable  :: error_diis(:,:)
+  double precision,allocatable  :: F_diis(:,:)
+  double precision,allocatable  :: c(:,:)
+  double precision,allocatable  :: cp(:,:)
+  double precision,allocatable  :: eGT(:)
+  double precision,allocatable  :: eOld(:)
+  double precision,allocatable  :: Omega1s(:),Omega1t(:)
+  double precision,allocatable  :: X1s(:,:),X1t(:,:)
+  double precision,allocatable  :: Y1s(:,:),Y1t(:,:)
+  double precision,allocatable  :: rho1s(:,:,:),rho1t(:,:,:)
+  double precision,allocatable  :: Omega2s(:),Omega2t(:)
+  double precision,allocatable  :: X2s(:,:),X2t(:,:)
+  double precision,allocatable  :: Y2s(:,:),Y2t(:,:)
+  double precision,allocatable  :: rho2s(:,:,:),rho2t(:,:,:)
+  double precision,allocatable  :: P(:,:)
+  double precision,allocatable  :: F(:,:)
+  double precision,allocatable  :: Fp(:,:)
+  double precision,allocatable  :: J(:,:)
+  double precision,allocatable  :: K(:,:)
+  double precision,allocatable  :: SigT(:,:)
+  double precision,allocatable  :: SigTp(:,:)
+  double precision,allocatable  :: SigTm(:,:)
+  double precision,allocatable  :: Z(:)
+  double precision,allocatable  :: error(:,:)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|       Self-consistent qsGT calculation       |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Dimensions of the pp-RPA linear reponse matrices
+
+  nOOs = nO*nO
+  nVVs = nV*nV
+
+  nOOt = nO*(nO - 1)/2
+  nVVt = nV*(nV - 1)/2
+
+! Warning 
+
+  write(*,*) '!! ERIs in MO basis will be overwritten in qsGT !!'
+  write(*,*)
+
+! Stuff 
+
+  nBasSq = nBas*nBas
+
+! TDA for T
+
+  if(TDA_T) then 
+    write(*,*) 'Tamm-Dancoff approximation for T-matrix!'
+    write(*,*)
+  end if
+
+! TDA 
+
+  if(TDA) then 
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*)
+  end if
+
+! Memory allocation
+
+  allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
+           J(nBas,nBas),K(nBas,nBas),SigT(nBas,nBas),SigTp(nBas,nBas),SigTm(nBas,nBas),Z(nBas),     & 
+           error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
+
+  allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),        &
+           Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs),        &
+           rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs),        &
+           Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),        &
+           Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt),        &
+           rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
+
+! Initialization
+  
+  nSCF            = -1
+  n_diis          = 0
+  ispin           = 1
+  Conv            = 1d0
+  P(:,:)          = PHF(:,:)
+  eGT(:)          = eHF(:)
+  eOld(:)         = eHF(:)
+  c(:,:)          = cHF(:,:)
+  F_diis(:,:)     = 0d0
+  error_diis(:,:) = 0d0
+  rcond           = 1d0
+
+!------------------------------------------------------------------------
+! Main loop
+!------------------------------------------------------------------------
+
+  do while(Conv > thresh .and. nSCF <= maxSCF)
+
+    ! Increment
+
+    nSCF = nSCF + 1
+
+    ! Buid Coulomb matrix
+
+    call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J)
+
+    ! Compute exchange part of the self-energy 
+
+    call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
+
+    ! AO to MO transformation of two-electron integrals
+
+    call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
+
+    ! Compute linear response
+
+    ispin  = 1
+    iblock = 3
+
+    call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO,  &
+                            Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
+
+    ispin  = 2
+    iblock = 4
+
+    call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO,  &
+                          Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
+
+    ! Compute correlation part of the self-energy 
+
+    SigT(:,:) = 0d0
+    Z(:)      = 0d0
+
+    iblock =  3
+
+    call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, &
+                                    X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+
+    call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
+                             Omega1s,rho1s,Omega2s,rho2s,SigT)
+
+    call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
+                                        Omega1s,rho1s,Omega2s,rho2s,Z)
+
+    iblock =  4
+
+    call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, &
+                                    X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+
+    call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
+                             Omega1t,rho1t,Omega2t,rho2t,SigT)
+
+    call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
+                                        Omega1t,rho1t,Omega2t,rho2t,Z)
+
+    Z(:) = 1d0/(1d0 - Z(:))
+
+    ! Make correlation self-energy Hermitian and transform it back to AO basis
+   
+    SigTp = 0.5d0*(SigT + transpose(SigT))
+    SigTm = 0.5d0*(SigT - transpose(SigT))
+
+    call MOtoAO_transform(nBas,S,c,SigTp)
+ 
+    ! Solve the quasi-particle equation
+
+    F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigTp(:,:)
+
+    ! Compute commutator and convergence criteria
+
+    error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
+
+    ! DIIS extrapolation 
+
+    n_diis = min(n_diis+1,max_diis)
+    if(abs(rcond) > 1d-7) then 
+      call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
+    else
+      n_diis = 0
+    end if
+
+    ! Diagonalize Hamiltonian in AO basis
+
+    Fp = matmul(transpose(X),matmul(F,X))
+    cp(:,:) = Fp(:,:)
+    call diagonalize_matrix(nBas,cp,eGT)
+    c = matmul(X,cp)
+    SigTp = matmul(transpose(c),matmul(SigTp,c))
+
+    ! Compute new density matrix in the AO basis
+
+    P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
+
+    ! Save quasiparticles energy for next cycle
+
+    Conv = maxval(abs(eGT - eOld))
+    eOld(:) = eGT(:)
+
+    !------------------------------------------------------------------------
+    !   Compute total energy
+    !------------------------------------------------------------------------
+
+    ! Kinetic energy
+
+    ET = trace_matrix(nBas,matmul(P,T))
+
+    ! Potential energy
+
+    EV = trace_matrix(nBas,matmul(P,V))
+
+    ! Coulomb energy
+
+    EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
+
+    ! Exchange energy
+
+    Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
+
+    ! Total energy
+
+    EqsGT = ET + EV + EJ + Ex 
+
+    ! Print results
+
+    call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
+    call print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigTp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
+
+  enddo
+!------------------------------------------------------------------------
+! End main loop
+!------------------------------------------------------------------------
+
+! Compute the ppRPA correlation energy
+
+  ispin  = 1
+  iblock = 3
+  call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO,  &
+                          Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
+  ispin  = 2
+  iblock = 4
+  call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO,  &
+                          Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
+  EcRPA(1) = EcRPA(1) - EcRPA(2)
+  EcRPA(2) = 3d0*EcRPA(2)
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (singlet) =',EcRPA(1)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (triplet) =',EcRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy           =',EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
+! Did it actually converge?
+
+  if(nSCF == maxSCF+1) then
+
+    write(*,*)
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)'                 Convergence failed                 '
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)
+
+    stop
+
+  endif
+
+! Deallocate memory
+
+  deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis)
+
+! Perform BSE calculation
+
+  if(BSE) then
+
+    call Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,   &
+                                Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, &
+                                ERI_MO,dipole_int_MO,eGT,eGT,EcBSE)
+
+    if(exchange_kernel) then
+
+      EcBSE(1) = 0.5d0*EcBSE(1)
+      EcBSE(2) = 1.5d0*EcBSE(2)
+
+    end if
+
+    write(*,*)
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT correlation energy (singlet) =',EcBSE(1)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT correlation energy (triplet) =',EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT correlation energy           =',EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT total energy                 =',ENuc + EqsGT + EcBSE(1) + EcBSE(2)
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,*)
+
+!   Compute the BSE correlation energy via the adiabatic connection 
+
+    if(doACFDT) then
+
+      write(*,*) '------------------------------------------------------'
+      write(*,*) 'Adiabatic connection version of BSE correlation energy'
+      write(*,*) '------------------------------------------------------'
+      write(*,*)
+
+      if(doXBS) then
+
+        write(*,*) '*** scaled screening version (XBS) ***'
+        write(*,*)
+
+      end if
+
+      call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC)
+
+      write(*,*)
+      write(*,*)'-------------------------------------------------------------------------------'
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT correlation energy (singlet) =',EcAC(1)
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT correlation energy (triplet) =',EcAC(2)
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT correlation energy           =',EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT total energy                 =',ENuc + EqsGT + EcAC(1) + EcAC(2)
+      write(*,*)'-------------------------------------------------------------------------------'
+      write(*,*)
+
+    end if
+
+  end if
+
+end subroutine qsGT
diff --git a/src/MBPT/self_energy_Tmatrix.f90 b/src/MBPT/self_energy_Tmatrix.f90
index 651b101..fc06791 100644
--- a/src/MBPT/self_energy_Tmatrix.f90
+++ b/src/MBPT/self_energy_Tmatrix.f90
@@ -28,7 +28,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
 
 ! Output variables
 
-  double precision,intent(out)  :: SigT(nBas)
+  double precision,intent(inout)  :: SigT(nBas,nBas)
 
 !----------------------------------------------
 ! Occupied part of the T-matrix self-energy 
@@ -38,8 +38,8 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
     do q=nC+1,nBas-nR
       do i=nC+1,nO
         do cd=1,nVV
-          eps     = e(p)    + e(i) - Omega1(cd)
-          SigT(p) = SigT(p) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
+          eps       = e(p)      + e(i) - Omega1(cd)
+          SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
         enddo
       enddo
     enddo
@@ -51,10 +51,10 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
 
   do p=nC+1,nBas-nR
     do q=nC+1,nBas-nR
-      do a=1,nV-nR
+      do a=nO+1,nBas-nR
         do kl=1,nOO
-          eps     = e(p)    + e(nO+a) - Omega2(kl)
-          SigT(p) = SigT(p) + rho2(p,nO+a,kl)*rho2(q,nO+a,kl)*eps/(eps**2 + eta**2)
+          eps       = e(p)      + e(a) - Omega2(kl)
+          SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2)
         enddo
       enddo
     enddo
diff --git a/src/MBPT/self_energy_Tmatrix_diag.f90 b/src/MBPT/self_energy_Tmatrix_diag.f90
index 793c52e..f18ac6e 100644
--- a/src/MBPT/self_energy_Tmatrix_diag.f90
+++ b/src/MBPT/self_energy_Tmatrix_diag.f90
@@ -28,7 +28,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
 
 ! Output variables
 
-  double precision,intent(out)  :: SigT(nBas)
+  double precision,intent(inout)  :: SigT(nBas)
 
 !----------------------------------------------
 ! Occupied part of the T-matrix self-energy 
@@ -48,10 +48,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
 !----------------------------------------------
 
   do p=nC+1,nBas-nR
-    do a=1,nV-nR
+    do a=nO+1,nBas-nR
       do kl=1,nOO
-        eps     = e(p)    + e(nO+a) - Omega2(kl)
-        SigT(p) = SigT(p) + rho2(p,nO+a,kl)**2*eps/(eps**2 + eta**2)
+        eps     = e(p)    + e(a) - Omega2(kl)
+        SigT(p) = SigT(p) + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 5ed4a55..dbf2dd0 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -1062,6 +1062,26 @@ program QuAcK
 
   end if
 
+!------------------------------------------------------------------------
+! Perform qsGT calculation
+!------------------------------------------------------------------------
+
+  if(doqsGT) then 
+
+    call cpu_time(start_qsGT)
+
+    call qsGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,         &
+              BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc, & 
+              nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+
+    call cpu_time(end_qsGT)
+
+    t_qsGT = end_qsGT - start_qsGT
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_qsGT,' seconds'
+    write(*,*)
+
+  end if
+
 !------------------------------------------------------------------------
 ! Information for Monte Carlo calculations
 !------------------------------------------------------------------------