diff --git a/input/methods b/input/methods
index 5921001..6e62625 100644
--- a/input/methods
+++ b/input/methods
@@ -1,5 +1,5 @@
 # RHF UHF GHF ROHF 
-  T   F   T   F   
+  T   F   F   F   
 # MP2 MP3 
   T   T   
 # CCD pCCD DCD CCSD CCSD(T) 
@@ -13,8 +13,8 @@
 # G0F2 evGF2 qsGF2 G0F3 evGF3
   T    T     F     F    F
 # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
-  F    F    F    F        F      F
+  T    T    T    F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
   F      F      F      F      F      F
 # Rtest Utest Gtest
-  T     F     T
+  T     F     F
diff --git a/input/options b/input/options
index 980f0d4..92084cd 100644
--- a/input/options
+++ b/input/options
@@ -9,9 +9,9 @@
 # GF: maxSCF thresh  DIIS lin eta renorm reg
       256    0.00001 5    F   0.0 0      F
 # GW: maxSCF thresh  DIIS lin eta TDA_W reg
-      10   0.00001 5    F   0.0 F     F 
+      256   0.00001  5    F   0.0 F     F 
 # GT: maxSCF thresh  DIIS lin eta TDA_T reg
-      256    0.00001 5     F  0.0 F     F  
+      256    0.00001 5    F   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: phBSE phBSE2 ppBSE dBSE dTDA
diff --git a/src/GT/G0T0eh.f90 b/src/GT/RG0T0eh.f90
similarity index 88%
rename from src/GT/G0T0eh.f90
rename to src/GT/RG0T0eh.f90
index 055144e..f44c506 100644
--- a/src/GT/G0T0eh.f90
+++ b/src/GT/RG0T0eh.f90
@@ -1,6 +1,5 @@
-subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & 
-                  singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,     & 
-                  ERI,dipole_int,eHF)
+subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & 
+                   singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform ehG0T0 calculation
 
@@ -10,6 +9,8 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
   logical,intent(in)            :: doXBS
@@ -73,9 +74,9 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|        One-shot G0T0eh calculation           |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Restricted G0T0eh Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! Initialization
@@ -167,4 +168,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
 
   call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','RG0T0eh correlation energy',EcRPA)
+    call dump_test_value('R','RG0T0eh HOMO energy',eGT(nO))
+    call dump_test_value('R','RG0T0eh LUMO energy',eGT(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GT/G0T0pp.f90 b/src/GT/RG0T0pp.f90
similarity index 92%
rename from src/GT/G0T0pp.f90
rename to src/GT/RG0T0pp.f90
index 7a1e94c..9e21b68 100644
--- a/src/GT/G0T0pp.f90
+++ b/src/GT/RG0T0pp.f90
@@ -1,5 +1,5 @@
-subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 
-                  linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 
+                   linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
 
@@ -8,6 +8,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
   logical,intent(in)            :: doXBS
@@ -65,9 +67,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|          One-shot G0T0pp calculation          |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Restricted G0T0pp Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 
@@ -260,8 +262,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy           =',EcBSE(1) + EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy           =',sum(EcBSE),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp total energy                 =',ENuc + ERHF + sum(EcBSE),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -296,8 +298,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
       write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
-      write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy           =',EcBSE(1) + EcBSE(2),' au'
-      write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
+      write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy           =',sum(EcBSE),' au'
+      write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp total energy                 =',ENuc + ERHF + sum(EcBSE),' au'
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,*)
 
@@ -315,11 +317,21 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy =',EcBSE(1) + EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy =',sum(EcBSE),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+! Testing zone
+
+  if(dotest) then
+  
+    call dump_test_value('R','RG0T0pp correlation energy',sum(EcRPA))
+    call dump_test_value('R','RG0T0pp HOMO energy',eGT(nO))
+    call dump_test_value('R','RG0T0pp LUMO energy',eGT(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90
index 22eddf0..e35890c 100644
--- a/src/GT/RGT.f90
+++ b/src/GT/RGT.f90
@@ -1,4 +1,4 @@
-subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
+subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
                exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,   &
                linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,      &
                ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
@@ -10,12 +10,14 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
 
 ! Input variables
 
-  logical                       :: doG0T0pp
-  logical                       :: doevGTpp
-  logical                       :: doqsGTpp
-  logical                       :: doG0T0eh
-  logical                       :: doevGTeh
-  logical                       :: doqsGTeh
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doG0T0pp
+  logical,intent(in)            :: doevGTpp
+  logical,intent(in)            :: doqsGTpp
+  logical,intent(in)            :: doG0T0eh
+  logical,intent(in)            :: doevGTeh
+  logical,intent(in)            :: doqsGTeh
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -74,7 +76,7 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doG0T0pp) then
     
     call wall_time(start_GT)
-    call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
+    call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GT)
   
@@ -91,8 +93,8 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doevGTpp) then
     
     call wall_time(start_GT)
-    call evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
-                linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
+                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GT)
   
     t_GT = end_GT - start_GT
@@ -108,9 +110,9 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doqsGTpp) then 
 
     call wall_time(start_GT)
-    call qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,          &
-                eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, & 
-                PHF,cHF,epsHF)
+    call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,          &
+                 eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, & 
+                 PHF,cHF,epsHF)
     call wall_time(end_GT)
 
     t_GT = end_GT - start_GT
@@ -126,8 +128,8 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doG0T0eh) then
     
     call wall_time(start_GT)
-    call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
-                linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
+                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GT)
   
     t_GT = end_GT - start_GT
@@ -143,8 +145,8 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doevGTeh) then
 
     call wall_time(start_GT)
-    call evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
-                singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
+                 singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GT)
 
     t_GT = end_GT - start_GT
@@ -160,9 +162,9 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doqsGTeh) then 
 
     call wall_time(start_GT)
-    call qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
-                eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, & 
-                PHF,cHF,epsHF)
+    call qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
+                 eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, & 
+                 PHF,cHF,epsHF)
     call wall_time(end_GT)
 
     t_GT = end_GT - start_GT
diff --git a/src/GT/UG0T0pp.f90 b/src/GT/UG0T0pp.f90
index a34d926..99de73a 100644
--- a/src/GT/UG0T0pp.f90
+++ b/src/GT/UG0T0pp.f90
@@ -1,4 +1,4 @@
-subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
+subroutine UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
                    spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, &
                    nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
                    dipole_int_aa,dipole_int_bb,cHF,eHF)
@@ -10,6 +10,8 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
   logical,intent(in)            :: doXBS
@@ -60,17 +62,16 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
   double precision,allocatable  :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:)
   double precision,allocatable  :: SigT(:,:)
   double precision,allocatable  :: Z(:,:)
-  double precision,allocatable  :: eG0T0(:,:)
+  double precision,allocatable  :: eGT(:,:)
 
 ! Output variables
 
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|          One-shot G0T0 calculation           |'
-  write(*,*)'|         *** Unrestricted version ***         |'
-  write(*,*)'************************************************'
+  write(*,*)'***********************************'
+  write(*,*)'* Unrestricted G0T0pp Calculation *'
+  write(*,*)'***********************************'
   write(*,*)
 
 ! Dimensions of the pp-URPA linear reponse matrices
@@ -101,7 +102,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
            Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb),   &
            Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb),   &
            rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), &
-           SigT(nBas,nspin),Z(nBas,nspin),eG0T0(nBas,nspin))
+           SigT(nBas,nspin),Z(nBas,nspin),eGT(nBas,nspin))
 
 !----------------------------------------------
 ! alpha-beta block
@@ -183,7 +184,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
 
   if(linearize) then 
 
-    eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:)
+    eGT(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:)
 
   else
  
@@ -203,7 +204,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
   ispin  = 1
   iblock = 3 
 
-  call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, &
+  call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
 
 !alpha-alpha block
@@ -211,7 +212,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
   ispin  = 2
   iblock = 4 
   
-  call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, &
+  call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) 
  
   Ecaa = EcRPA(2)
@@ -220,7 +221,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
  
   iblock = 7
   
-  call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, &
+  call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
   
   Ecbb = EcRPA(2) 
@@ -228,7 +229,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
   EcRPA(1) = EcRPA(1) - EcRPA(2)
   EcRPA(2) = 3d0*EcRPA(2)
 
-  call print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eG0T0,EcGM,EcRPA)
+  call print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
 
 ! Free memory
 
@@ -236,4 +237,16 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
              Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
              Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb)
 
+! Testing zone
+
+  if(dotest) then
+  
+    call dump_test_value('U','UG0T0pp correlation energy',sum(EcRPA))
+    call dump_test_value('U','UG0T0pp HOMOa energy',eGT(nO(1),1))
+    call dump_test_value('U','UG0T0pp LUMOa energy',eGT(nO(1)+1,1))
+    call dump_test_value('U','UG0T0pp HOMOa energy',eGT(nO(2),2))
+    call dump_test_value('U','UG0T0pp LUMOa energy',eGT(nO(2)+1,2))
+
+  end if
+
 end subroutine 
diff --git a/src/GT/UGT.f90 b/src/GT/UGT.f90
index 15c43ea..2458346 100644
--- a/src/GT/UGT.f90
+++ b/src/GT/UGT.f90
@@ -1,4 +1,4 @@
-subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT,        &
+subroutine UGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
                exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
                linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,             &
                ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
@@ -10,12 +10,14 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
 
 ! Input variables
 
-  logical                       :: doG0T0pp
-  logical                       :: doevGTpp
-  logical                       :: doqsGTpp
-  logical                       :: doG0T0eh
-  logical                       :: doevGTeh
-  logical                       :: doqsGTeh
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doG0T0pp
+  logical,intent(in)            :: doevGTpp
+  logical,intent(in)            :: doqsGTpp
+  logical,intent(in)            :: doG0T0eh
+  logical,intent(in)            :: doevGTeh
+  logical,intent(in)            :: doqsGTeh
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -77,7 +79,7 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doG0T0pp) then
     
     call wall_time(start_GT)
-    call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,      &
+    call UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,      &
                  linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
                  dipole_int_aa,dipole_int_bb,cHF,epsHF)
     call wall_time(end_GT)
@@ -95,7 +97,7 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doevGTpp) then
     
     call wall_time(start_GT)
-    call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,    &
+    call evUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
                  linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & 
                  cHF,epsHF)
     call wall_time(end_GT)
@@ -113,7 +115,7 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
   if(doqsGTpp) then 
 
     call wall_time(start_GT)
-    call qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & 
+    call qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & 
                  eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,   &
                  dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) 
     call wall_time(end_GT)
diff --git a/src/GT/evGTeh.f90 b/src/GT/evRGTeh.f90
similarity index 93%
rename from src/GT/evGTeh.f90
rename to src/GT/evRGTeh.f90
index b316755..b8b8e3c 100644
--- a/src/GT/evGTeh.f90
+++ b/src/GT/evRGTeh.f90
@@ -1,5 +1,5 @@
-subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & 
-                  singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & 
+                   singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform self-consistent eigenvalue-only ehGT calculation
 
@@ -8,6 +8,8 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -70,9 +72,9 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|    Self-consistent evGTeh calculation        |'
-  write(*,*)'************************************************'
+  write(*,*)'**********************************'
+  write(*,*)'* Restricted evRGTeh Calculation *'
+  write(*,*)'**********************************'
   write(*,*)
 
 ! TDA for T
@@ -296,4 +298,14 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
 ! end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','evRGTeh correlation energy',EcRPA)
+    call dump_test_value('R','evRGTeh HOMO energy',eGT(nO))
+    call dump_test_value('R','evRGTeh LUMO energy',eGT(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GT/evGTpp.f90 b/src/GT/evRGTpp.f90
similarity index 92%
rename from src/GT/evGTpp.f90
rename to src/GT/evRGTpp.f90
index c33bf4b..97ef5dd 100644
--- a/src/GT/evGTpp.f90
+++ b/src/GT/evRGTpp.f90
@@ -1,5 +1,5 @@
-subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & 
-                  linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & 
+                   linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
 
@@ -8,6 +8,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -74,9 +76,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|      Self-consistent evGTpp calculation      |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Restricted evGTpp Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! TDA for T
@@ -265,8 +267,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1)
     write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy           =',sum(EcBSE)
+    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp total energy                 =',ENuc + ERHF + sum(EcBSE)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -301,8 +303,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1)
       write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (triplet) =',EcBSE(2)
-      write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy           =',EcBSE(1) + EcBSE(2)
-      write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
+      write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy           =',sum(EcBSE)
+      write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp total energy                 =',ENuc + ERHF + sum(EcBSE)
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,*)
 
@@ -310,4 +312,15 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
 
   end if
 
+
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','evRGTpp correlation energy',sum(EcRPA))
+    call dump_test_value('R','evRGTpp HOMO energy',eGT(nO))
+    call dump_test_value('R','evRGTpp LUMO energy',eGT(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GT/evUGTpp.f90 b/src/GT/evUGTpp.f90
index 5217ace..f69bc96 100644
--- a/src/GT/evUGTpp.f90
+++ b/src/GT/evUGTpp.f90
@@ -1,4 +1,4 @@
-subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
+subroutine evUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
                    TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
                    eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa, &
                    ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF)
@@ -9,6 +9,9 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
   include 'parameters.h'
 
 ! Input variables
+
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -73,9 +76,9 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|      Self-consistent evUGT calculation        |'
-  write(*,*)'************************************************'
+  write(*,*)'***********************************'
+  write(*,*)'* Unrestricted evGTpp Calculation *'
+  write(*,*)'***********************************'
   write(*,*)
 
 ! Dimensions of the pp-URPA linear reponse matrices
@@ -284,4 +287,16 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
                Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
                Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb)
 
+! Testing zone
+  
+  if(dotest) then
+  
+    call dump_test_value('U','evUGTpp correlation energy',sum(EcRPA))
+    call dump_test_value('U','evUGTpp HOMOa energy',eGT(nO(1),1))
+    call dump_test_value('U','evUGTpp LUMOa energy',eGT(nO(1)+1,1))
+    call dump_test_value('U','evUGTpp HOMOa energy',eGT(nO(2),2))
+    call dump_test_value('U','evUGTpp LUMOa energy',eGT(nO(2)+1,2))
+  
+  end if
+
 end subroutine 
diff --git a/src/GT/qsGTeh.f90 b/src/GT/qsRGTeh.f90
similarity index 87%
rename from src/GT/qsGTeh.f90
rename to src/GT/qsRGTeh.f90
index 7b7ff6c..9920489 100644
--- a/src/GT/qsGTeh.f90
+++ b/src/GT/qsRGTeh.f90
@@ -1,6 +1,6 @@
-subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,       & 
-                  dBSE,dTDA,singlet,triplet,eta,regularize,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)
+subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA, & 
+                   dBSE,dTDA,singlet,triplet,eta,regularize,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 GTeh calculation
 
@@ -9,6 +9,8 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -73,7 +75,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
   double precision              :: dipole(ncart)
 
   logical                       :: print_T = .true.
-  double precision,allocatable  :: error_diis(:,:)
+  double precision,allocatable  :: err_diis(:,:)
   double precision,allocatable  :: F_diis(:,:)
   double precision,allocatable  :: Aph(:,:)
   double precision,allocatable  :: Bph(:,:)
@@ -94,14 +96,14 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
   double precision,allocatable  :: Sig(:,:)
   double precision,allocatable  :: Sigp(:,:)
   double precision,allocatable  :: Z(:)
-  double precision,allocatable  :: error(:,:)
+  double precision,allocatable  :: err(:,:)
 
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|     Self-consistent qsGTeh calculation       |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Restricted qsGTeh Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! Warning 
@@ -131,21 +133,21 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
   allocate(Aph(nS,nS),Bph(nS,nS),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),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & 
-           rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
+           rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
 
 ! Initialization
   
-  nSCF            = -1
-  n_diis          = 0
-  ispin           = 2
-  Conv            = 1d0
-  P(:,:)          = PHF(:,:)
-  eGT(:)          = eHF(:)
-  eOld(:)         = eHF(:)
-  c(:,:)          = cHF(:,:)
-  F_diis(:,:)     = 0d0
-  error_diis(:,:) = 0d0
-  rcond           = 0d0
+  nSCF          = -1
+  n_diis        = 0
+  ispin         = 2
+  Conv          = 1d0
+  P(:,:)        = PHF(:,:)
+  eGT(:)        = eHF(:)
+  eOld(:)       = eHF(:)
+  c(:,:)        = cHF(:,:)
+  F_diis(:,:)   = 0d0
+  err_diis(:,:) = 0d0
+  rcond         = 0d0
 
 !------------------------------------------------------------------------
 ! Main loop
@@ -198,14 +200,14 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
     ! Compute commutator and convergence criteria
 
-    error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
+    err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
 
     ! DIIS extrapolation 
 
     if(max_diis > 1) then
 
       n_diis = min(n_diis+1,max_diis)
-      call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
+      call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
 
     end if
 
@@ -223,7 +225,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
     ! Save quasiparticles energy for next cycle
 
-    Conv = maxval(abs(error))
+    Conv = maxval(abs(err))
     eOld(:) = eGT(:)
 
     !------------------------------------------------------------------------
@@ -276,7 +278,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
 ! Deallocate memory
 
-  deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,error,error_diis,F_diis)
+  deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,err,err_diis,F_diis)
 
 ! Perform BSE calculation
 
@@ -332,4 +334,14 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
 ! end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','qsRGTeh correlation energy',EcRPA)
+    call dump_test_value('R','qsRGTeh HOMO energy',eGT(nO))
+    call dump_test_value('R','qsRGTeh LUMO energy',eGT(nO+1))
+
+  end if
+
 end subroutine
diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsRGTpp.f90
similarity index 93%
rename from src/GT/qsGTpp.f90
rename to src/GT/qsRGTpp.f90
index 118ae13..619fc03 100644
--- a/src/GT/qsGTpp.f90
+++ b/src/GT/qsRGTpp.f90
@@ -1,6 +1,6 @@
-subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,          & 
-                  dBSE,dTDA,singlet,triplet,eta,regularize,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)
+subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,          & 
+                   dBSE,dTDA,singlet,triplet,eta,regularize,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
 
@@ -9,6 +9,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -98,9 +99,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|       Self-consistent qsGTpp calculation     |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Restricted qsGTpp Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! Dimensions of the pp-RPA linear reponse matrices
@@ -341,8 +342,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1)
     write(*,'(2X,A50,F20.10)') 'Tr@BphSE@qsGTpp correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp total energy                 =',ENuc + EqsGT + EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy           =',sum(EcBSE)
+    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp total energy                 =',ENuc + EqsGT + sum(EcBSE)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -370,8 +371,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1)
       write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy (triplet) =',EcBSE(2)
-      write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy           =',EcBSE(1) + EcBSE(2)
-      write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp total energy                 =',ENuc + EqsGT + EcBSE(1) + EcBSE(2)
+      write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy           =',sum(EcBSE)
+      write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp total energy                 =',ENuc + EqsGT + sum(EcBSE)
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,*)
 
@@ -379,4 +380,15 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
 
   end if
 
+
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','qsRGTpp correlation energy',sum(EcRPA))
+    call dump_test_value('R','qsRGTpp HOMO energy',eGT(nO))
+    call dump_test_value('R','qsRGTpp LUMO energy',eGT(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90
index 390f1dc..821f93f 100644
--- a/src/GT/qsUGTpp.f90
+++ b/src/GT/qsUGTpp.f90
@@ -1,4 +1,4 @@
-subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
+subroutine qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
                    TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
                    eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,&
                    ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
@@ -9,6 +9,9 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
   include 'parameters.h'
 
 ! Input variables
+
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -103,9 +106,9 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|      Self-consistent qsUGT calculation        |'
-  write(*,*)'************************************************'
+  write(*,*)'***********************************'
+  write(*,*)'* Unrestricted evGTpp Calculation *'
+  write(*,*)'***********************************'
   write(*,*)
 
 ! Dimensions of the pp-URPA linear reponse matrices
@@ -409,4 +412,16 @@ write(*,*) 'EcGM', EcGM(1)
 
   deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,Z,error,error_diis,F_diis) 
 
+! Testing zone
+  
+  if(dotest) then
+  
+    call dump_test_value('U','qsUGTpp correlation energy',sum(EcRPA))
+    call dump_test_value('U','qsUGTpp HOMOa energy',eGT(nO(1),1))
+    call dump_test_value('U','qsUGTpp LUMOa energy',eGT(nO(1)+1,1))
+    call dump_test_value('U','qsUGTpp HOMOa energy',eGT(nO(2),2))
+    call dump_test_value('U','qsUGTpp LUMOa energy',eGT(nO(2)+1,2))
+  
+  end if
+
 end subroutine 
diff --git a/test/methods.test b/test/methods.test
index 5921001..6e62625 100644
--- a/test/methods.test
+++ b/test/methods.test
@@ -1,5 +1,5 @@
 # RHF UHF GHF ROHF 
-  T   F   T   F   
+  T   F   F   F   
 # MP2 MP3 
   T   T   
 # CCD pCCD DCD CCSD CCSD(T) 
@@ -13,8 +13,8 @@
 # G0F2 evGF2 qsGF2 G0F3 evGF3
   T    T     F     F    F
 # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
-  F    F    F    F        F      F
+  T    T    T    F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
   F      F      F      F      F      F
 # Rtest Utest Gtest
-  T     F     T
+  T     F     F
diff --git a/test/options.test b/test/options.test
index 980f0d4..92084cd 100644
--- a/test/options.test
+++ b/test/options.test
@@ -9,9 +9,9 @@
 # GF: maxSCF thresh  DIIS lin eta renorm reg
       256    0.00001 5    F   0.0 0      F
 # GW: maxSCF thresh  DIIS lin eta TDA_W reg
-      10   0.00001 5    F   0.0 F     F 
+      256   0.00001  5    F   0.0 F     F 
 # GT: maxSCF thresh  DIIS lin eta TDA_T reg
-      256    0.00001 5     F  0.0 F     F  
+      256    0.00001 5    F   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: phBSE phBSE2 ppBSE dBSE dTDA