diff --git a/input/methods b/input/methods
index 682b927..5921001 100644
--- a/input/methods
+++ b/input/methods
@@ -1,20 +1,20 @@
 # RHF UHF GHF ROHF 
-  T   T   T   T   
+  T   F   T   F   
 # MP2 MP3 
   T   T   
 # CCD pCCD DCD CCSD CCSD(T) 
-  F   F    F   F    F
+  T   T    T   T    F
 # drCCD rCCD crCCD lCCD
-  F     F    F     F
+  T     T    T     T
 # CIS CIS(D) CID CISD FCI
-  F   F      F   F    F
+  T   F      F   F    F
 # phRPA phRPAx crRPA ppRPA 
-  T     T      F     T 
+  T     T      T     T 
 # G0F2 evGF2 qsGF2 G0F3 evGF3
-  F    F     F     F    F
+  T    T     F     F    F
 # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
   F    F    F    F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
   F      F      F      F      F      F
 # Rtest Utest Gtest
-  T     T     T
+  T     F     T
diff --git a/input/options b/input/options
index 979a4b1..980f0d4 100644
--- a/input/options
+++ b/input/options
@@ -4,8 +4,8 @@
       F
 # CC: maxSCF thresh     DIIS
       64     0.0000001  5
-# spin: TDA spin_conserved spin_flip 
-        F   T               T
+# spin: TDA singlet triplet
+        F   T       T
 # 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
diff --git a/src/CC/CCD.f90 b/src/CC/CCD.f90
index e7844ab..c44e00e 100644
--- a/src/CC/CCD.f90
+++ b/src/CC/CCD.f90
@@ -1,4 +1,4 @@
-subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine CCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! CCD module
 
@@ -6,6 +6,8 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -29,7 +31,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
   integer                       :: nSCF
   double precision              :: Conv
   double precision              :: EcMP2,EcMP3,EcMP4
-  double precision              :: ECCD,EcCCD
+  double precision              :: ECC,EcCC
   double precision,allocatable  :: seHF(:)
   double precision,allocatable  :: sERI(:,:,:,:)
   double precision,allocatable  :: dbERI(:,:,:,:)
@@ -191,13 +193,13 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
 
 !   Compute correlation energy
 
-    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCCD)
+    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCC)
 
     if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t,v,delta_OOVV,EcMP3)
 
 !   Dump results
 
-    ECCD = ERHF + EcCCD
+    ECC = ERHF + EcCC
 
     ! DIIS extrapolation
 
@@ -209,7 +211,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
     if(abs(rcond) < 1d-15) n_diis = 0
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
   enddo
   write(*,*)'----------------------------------------------------'
@@ -235,8 +237,8 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
   write(*,*)'----------------------------------------------------'
   write(*,*)'              CCD energy                            '
   write(*,*)'----------------------------------------------------'
-  write(*,'(1X,A30,1X,F15.10)')' E(CCD) = ',ECCD
-  write(*,'(1X,A30,1X,F15.10)')' Ec(CCD) = ',EcCCD
+  write(*,'(1X,A30,1X,F15.10)')' E(CCD) = ',ECC
+  write(*,'(1X,A30,1X,F15.10)')' Ec(CCD) = ',EcCC
   write(*,*)'----------------------------------------------------'
   write(*,*)
 
@@ -272,4 +274,12 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
 
   if(do_DIP_EOM_CC_2h) call DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t)
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','CCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/CCSD.f90 b/src/CC/CCSD.f90
index 42d0a90..0c895b4 100644
--- a/src/CC/CCSD.f90
+++ b/src/CC/CCSD.f90
@@ -1,4 +1,4 @@
-subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! CCSD module
 
@@ -6,6 +6,8 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,EN
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -302,4 +304,12 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,EN
 
   end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','CCSD correlation energy',EcCCSD)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/DCD.f90 b/src/CC/DCD.f90
index e7d856a..91344e8 100644
--- a/src/CC/DCD.f90
+++ b/src/CC/DCD.f90
@@ -1,4 +1,4 @@
-subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
+subroutine DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! DCD module
 
@@ -6,6 +6,8 @@ subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -24,7 +26,7 @@ subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
   integer                       :: nSCF
   double precision              :: Conv
   double precision              :: EcMP2,EcMP3,EcMP4
-  double precision              :: EDCD,EcDCD
+  double precision              :: ECC,EcCC
 
   double precision,allocatable  :: eO(:)
   double precision,allocatable  :: eV(:)
@@ -123,13 +125,13 @@ subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
     ! Compute correlation energy
 
-    EcDCD = 0d0
+    EcCC = 0d0
     do i=1,nO-nC
       do j=1,nO-nC
         do a=1,nV-nR
           do b=1,nV-nR
  
-            EcDCD = EcDCD + (2d0*OOVV(i,j,a,b) - OOVV(i,j,b,a))*t(i,j,a,b)
+            EcCC = EcCC + (2d0*OOVV(i,j,a,b) - OOVV(i,j,b,a))*t(i,j,a,b)
  
           enddo
         enddo
@@ -138,10 +140,10 @@ subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
     ! Dump results
 
-    EDCD = ERHF + EcDCD
+    ECC = ERHF + EcCC
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',EDCD+ENuc,'|',EcDCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
     ! Increment 
 
@@ -280,4 +282,12 @@ subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   endif
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','DCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90
index 2e5722d..5f2e531 100644
--- a/src/CC/RCC.f90
+++ b/src/CC/RCC.f90
@@ -1,4 +1,4 @@
-subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & 
+subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & 
                maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
 
 ! Coupled-cluster module
@@ -8,6 +8,8 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
 
 ! Input variables
 
+  logical                       :: dotest
+
   logical                       :: doCCD
   logical                       :: dopCCD
   logical                       :: doDCD
@@ -43,7 +45,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(doCCD) then
 
     call wall_time(start_CC)
-    call CCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -59,7 +61,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(doDCD) then
 
     call wall_time(start_CC)
-    call DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & 
+    call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & 
              ERI,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
@@ -78,7 +80,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(doCCSD) then
 
     call wall_time(start_CC)
-    call CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -94,7 +96,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(dodrCCD) then
 
     call wall_time(start_CC)
-    call drCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -110,7 +112,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(dorCCD) then
 
     call wall_time(start_CC)
-    call rCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,epsHF)
+    call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,epsHF)
     call wall_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -126,7 +128,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(docrCCD) then
 
     call wall_time(start_CC)
-    call crCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -142,7 +144,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(dolCCD) then
 
     call wall_time(start_CC)
-    call lCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -158,7 +160,7 @@ subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD,
   if(dopCCD) then
 
     call wall_time(start_CC)
-    call pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
 
     call wall_time(end_CC)
 
diff --git a/src/CC/crCCD.f90 b/src/CC/crCCD.f90
index d427ba7..ffc1881 100644
--- a/src/CC/crCCD.f90
+++ b/src/CC/crCCD.f90
@@ -1,4 +1,4 @@
-subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine crCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! Crossed-ring CCD module
 
@@ -6,6 +6,8 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -29,7 +31,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
   integer                       :: nSCF
   double precision              :: Conv
   double precision              :: EcMP2
-  double precision              :: ECCD,EcCCD
+  double precision              :: ECC,EcCC
   double precision,allocatable  :: seHF(:)
   double precision,allocatable  :: sERI(:,:,:,:)
   double precision,allocatable  :: dbERI(:,:,:,:)
@@ -155,11 +157,11 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
 
 !   Compute correlation energy
 
-    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
+    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC)
 
 !   Dump results
 
-    ECCD = ERHF + EcCCD
+    ECC = ERHF + EcCC
 
     ! DIIS extrapolation
 
@@ -171,7 +173,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
     if(abs(rcond) < 1d-15) n_diis = 0
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
   enddo
   write(*,*)'----------------------------------------------------'
@@ -197,9 +199,15 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
   write(*,*)'----------------------------------------------------'
   write(*,*)'       crossed-ring CCD energy                      '
   write(*,*)'----------------------------------------------------'
-  write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD  
-  write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD 
+  write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECC  
+  write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCC 
   write(*,*)'----------------------------------------------------'
   write(*,*)
 
+  if(dotest) then
+
+    call dump_test_value('R','crCCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/drCCD.f90 b/src/CC/drCCD.f90
index 33764d6..0a46948 100644
--- a/src/CC/drCCD.f90
+++ b/src/CC/drCCD.f90
@@ -1,4 +1,4 @@
-subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine drCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! Direct ring CCD module
 
@@ -6,6 +6,8 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -29,7 +31,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
   integer                       :: nSCF
   double precision              :: Conv
   double precision              :: EcMP2
-  double precision              :: ECCD,EcCCD
+  double precision              :: ECC,EcCC
   double precision,allocatable  :: seHF(:)
   double precision,allocatable  :: sERI(:,:,:,:)
 
@@ -147,12 +149,12 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
 
 !   Compute correlation energy
 
-    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
-    EcCCD = 2d0*EcCCD
+    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC)
+    EcCC = 2d0*EcCC
 
 !   Dump results
 
-    ECCD = ERHF + EcCCD
+    ECC = ERHF + EcCC
 
     ! DIIS extrapolation
 
@@ -164,7 +166,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
     if(abs(rcond) < 1d-15) n_diis = 0
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
   enddo
   write(*,*)'----------------------------------------------------'
@@ -190,9 +192,15 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
   write(*,*)'----------------------------------------------------'
   write(*,*)'              direct ring CCD energy                '
   write(*,*)'----------------------------------------------------'
-  write(*,'(1X,A30,1X,F15.10)')' E(drCCD) = ',ECCD  
-  write(*,'(1X,A30,1X,F15.10)')' Ec(drCCD) = ',EcCCD 
+  write(*,'(1X,A30,1X,F15.10)')' E(drCCD) = ',ECC  
+  write(*,'(1X,A30,1X,F15.10)')' Ec(drCCD) = ',EcCC 
   write(*,*)'----------------------------------------------------'
   write(*,*)
 
+  if(dotest) then
+
+    call dump_test_value('R','drCCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/lCCD.f90 b/src/CC/lCCD.f90
index b7c77ae..440933f 100644
--- a/src/CC/lCCD.f90
+++ b/src/CC/lCCD.f90
@@ -1,4 +1,4 @@
-subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine lCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! Ladder CCD module
 
@@ -6,6 +6,8 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -29,7 +31,7 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
   integer                       :: nSCF
   double precision              :: Conv
   double precision              :: EcMP2
-  double precision              :: ECCD,EcCCD
+  double precision              :: ECC,EcCC
   double precision,allocatable  :: seHF(:)
   double precision,allocatable  :: sERI(:,:,:,:)
   double precision,allocatable  :: dbERI(:,:,:,:)
@@ -168,11 +170,11 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
 
 !   Compute correlation energy
 
-    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
+    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC)
 
 !   Dump results
 
-    ECCD = ERHF + EcCCD
+    ECC = ERHF + EcCC
 
     ! DIIS extrapolation
 
@@ -184,7 +186,7 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
     if(abs(rcond) < 1d-15) n_diis = 0
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
   enddo
   write(*,*)'----------------------------------------------------'
@@ -210,9 +212,15 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
   write(*,*)'----------------------------------------------------'
   write(*,*)'              ladder CCD energy                     '
   write(*,*)'----------------------------------------------------'
-  write(*,'(1X,A30,1X,F15.10)')' E(lCCD)  = ',ECCD
-  write(*,'(1X,A30,1X,F15.10)')' Ec(lCCD) = ',EcCCD
+  write(*,'(1X,A30,1X,F15.10)')' E(lCCD)  = ',ECC
+  write(*,'(1X,A30,1X,F15.10)')' Ec(lCCD) = ',EcCC
   write(*,*)'----------------------------------------------------'
   write(*,*)
 
+  if(dotest) then
+
+    call dump_test_value('R','lCCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90
index 203db1a..916d9af 100644
--- a/src/CC/pCCD.f90
+++ b/src/CC/pCCD.f90
@@ -1,4 +1,4 @@
-subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
+subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! pair CCD module
 
@@ -6,6 +6,8 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -21,7 +23,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   integer                       :: nSCF
   double precision              :: Conv
-  double precision              :: ECCD,EcCCD
+  double precision              :: ECC,EcCC
 
   double precision,allocatable  :: eO(:)
   double precision,allocatable  :: eV(:)
@@ -91,10 +93,10 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   t(:,:) = -0.5d0*OOVV(:,:)/delta_OV(:,:)
 
-  EcCCD = 0d0
+  EcCC = 0d0
   do i=1,nO-nC
     do a=1,nV-nR
-      EcCCD = EcCCD + OOVV(i,a)*t(i,a)
+      EcCC = EcCC + OOVV(i,a)*t(i,a)
     end do
   end do
 
@@ -170,16 +172,16 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
    ! Compute correlation energy
 
-    EcCCD = 0d0
+    EcCC = 0d0
     do i=1,nO-nC    
       do a=1,nV-nR
-        EcCCD = EcCCD + OOVV(i,a)*t(i,a)
+        EcCC = EcCC + OOVV(i,a)*t(i,a)
      end do
    end do
 
    ! Dump results
 
-    ECCD = ERHF + EcCCD
+    ECC = ERHF + EcCC
 
     ! DIIS extrapolation
 
@@ -191,7 +193,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 !   if(abs(rcond) < 1d-15) n_diis = 0
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
   enddo
   write(*,*)'----------------------------------------------------'
@@ -212,4 +214,10 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   endif
 
+  if(dotest) then
+
+    call dump_test_value('R','pCCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90
index bba1607..ce61a66 100644
--- a/src/CC/rCCD.f90
+++ b/src/CC/rCCD.f90
@@ -1,4 +1,4 @@
-subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW)
+subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW)
 
 ! Ring CCD module
 
@@ -6,6 +6,8 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -31,7 +33,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
   integer                       :: nSCF
   double precision              :: Conv
   double precision              :: EcMP2
-  double precision              :: ECCD,EcCCD
+  double precision              :: ECC,EcCC
   double precision,allocatable  :: seHF(:)
   double precision,allocatable  :: seGW(:)
   double precision,allocatable  :: sERI(:,:,:,:)
@@ -161,11 +163,11 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
 
 !   Compute correlation energy
 
-    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCCD)
+    call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCC)
 
 !   Dump results
 
-    ECCD = ERHF + EcCCD
+    ECC = ERHF + EcCC
 
     ! DIIS extrapolation
 
@@ -177,7 +179,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
     if(abs(rcond) < 1d-15) n_diis = 0
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-      '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
 
   enddo
   write(*,*)'----------------------------------------------------'
@@ -203,8 +205,8 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
   write(*,*)'----------------------------------------------------'
   write(*,*)'              ring CCD energy                       '
   write(*,*)'----------------------------------------------------'
-  write(*,'(1X,A30,1X,F15.10)')' E(rCCD) = ',ECCD  
-  write(*,'(1X,A30,1X,F15.10)')' Ec(rCCD) = ',EcCCD 
+  write(*,'(1X,A30,1X,F15.10)')' E(rCCD) = ',ECC  
+  write(*,'(1X,A30,1X,F15.10)')' Ec(rCCD) = ',EcCC 
   write(*,*)'----------------------------------------------------'
   write(*,*)
 
@@ -223,4 +225,11 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
 
   if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
 
+
+  if(dotest) then
+
+    call dump_test_value('R','rCCD correlation energy',EcCC)
+
+  end if
+
 end subroutine 
diff --git a/src/CI/CID.f90 b/src/CI/CID.f90
index df14e0b..ef8d65e 100644
--- a/src/CI/CID.f90
+++ b/src/CI/CID.f90
@@ -1,4 +1,4 @@
-subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0)
+subroutine CID(dotest,singlet,triplet,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0)
 
 ! Perform configuration interaction with doubles
 
@@ -7,8 +7,10 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi
 
 ! Input variables
 
-  logical,intent(in)            :: singlet_manifold
-  logical,intent(in)            :: triplet_manifold
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: singlet
+  logical,intent(in)            :: triplet
   integer,intent(in)            :: nBasin
   integer,intent(in)            :: nCin
   integer,intent(in)            :: nOin
diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90
index fdd928e..38bce1d 100644
--- a/src/CI/CISD.f90
+++ b/src/CI/CISD.f90
@@ -1,4 +1,4 @@
-subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0)
+subroutine CISD(dotest,singlet,triplet,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0)
 
 ! Perform configuration interaction with singles and doubles
 
@@ -7,8 +7,10 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
 
 ! Input variables
 
-  logical,intent(in)            :: singlet_manifold
-  logical,intent(in)            :: triplet_manifold
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: singlet
+  logical,intent(in)            :: triplet
   integer,intent(in)            :: nBasin
   integer,intent(in)            :: nCin
   integer,intent(in)            :: nOin
diff --git a/src/CI/RCI.f90 b/src/CI/RCI.f90
index e6e21f7..84a7fe8 100644
--- a/src/CI/RCI.f90
+++ b/src/CI/RCI.f90
@@ -1,4 +1,4 @@
-subroutine RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, &
+subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, &
                epsHF,EHF,cHF,S)
 
 ! Configuration interaction module
@@ -8,11 +8,13 @@ subroutine RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR
 
 ! Input variables
 
-  logical                       :: doCIS
-  logical                       :: doCIS_D
-  logical                       :: doCID
-  logical                       :: doCISD
-  logical                       :: doFCI
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doCIS
+  logical,intent(in)            :: doCIS_D
+  logical,intent(in)            :: doCID
+  logical,intent(in)            :: doCISD
+  logical,intent(in)            :: doFCI
 
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
@@ -40,7 +42,7 @@ subroutine RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR
   if(doCIS) then
 
     call wall_time(start_CI)
-    call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
+    call RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
     call wall_time(end_CI)
 
     t_CI = end_CI - start_CI
@@ -56,7 +58,7 @@ subroutine RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR
   if(doCID) then
 
     call wall_time(start_CI)
-    call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
+    call CID(dotest,singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
     call wall_time(end_CI)
 
     t_CI = end_CI - start_CI
@@ -72,7 +74,7 @@ subroutine RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR
   if(doCISD) then
 
     call wall_time(start_CI)
-    call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
+    call CISD(dotest,singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
     call wall_time(end_CI)
 
     t_CI = end_CI - start_CI
diff --git a/src/CI/CIS.f90 b/src/CI/RCIS.f90
similarity index 87%
rename from src/CI/CIS.f90
rename to src/CI/RCIS.f90
index aff6b97..2d107a1 100644
--- a/src/CI/CIS.f90
+++ b/src/CI/RCIS.f90
@@ -1,4 +1,4 @@
-subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
+subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
 
 ! Perform configuration interaction single calculation`
 
@@ -7,6 +7,8 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   logical,intent(in)            :: doCIS_D
@@ -69,6 +71,14 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
     maxS = min(maxS,nS)
     if(doCIS_D) call CIS_D(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Om(1:maxS),A(:,1:maxS))
 
+    ! Testing zone
+  
+    if(dotest) then
+
+      call dump_test_value('R','RCIS singlet excitation energy',Om(1))
+
+    end if
+
   endif
 
   if(triplet) then
@@ -97,6 +107,14 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
     maxS = min(maxS,nS)
     if(doCIS_D) call CIS_D(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Om(1:maxS),A(:,1:maxS))
 
+    ! Testing zone
+ 
+    if(dotest) then
+ 
+      call dump_test_value('R','RCIS triplet excitation energy',Om(1))
+ 
+    end if
+
   endif
 
 end subroutine 
diff --git a/src/CI/UCI.f90 b/src/CI/UCI.f90
index 8ddae0e..ba2f3fb 100644
--- a/src/CI/UCI.f90
+++ b/src/CI/UCI.f90
@@ -1,4 +1,4 @@
-subroutine UCI(doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, & 
+subroutine UCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, & 
                ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,EHF,cHF,S,F)
 
 ! Configuration interaction module
@@ -8,11 +8,13 @@ subroutine UCI(doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC
 
 ! Input variables
 
-  logical                       :: doCIS
-  logical                       :: doCIS_D
-  logical                       :: doCID
-  logical                       :: doCISD
-  logical                       :: doFCI
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doCIS
+  logical,intent(in)            :: doCIS_D
+  logical,intent(in)            :: doCID
+  logical,intent(in)            :: doCISD
+  logical,intent(in)            :: doFCI
 
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
@@ -44,7 +46,7 @@ subroutine UCI(doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC
   if(doCIS) then
 
     call wall_time(start_CI)
-    call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb, & 
+    call UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb, & 
               ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
     call wall_time(end_CI)
 
diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90
index bd049ff..334c27c 100644
--- a/src/CI/UCIS.f90
+++ b/src/CI/UCIS.f90
@@ -1,4 +1,4 @@
-subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb, &
+subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb, &
                 dipole_int_aa,dipole_int_bb,eHF,cHF,S)
 
 ! Perform configuration interaction single calculation`
@@ -8,6 +8,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
   integer,intent(in)            :: nBas
@@ -88,6 +90,14 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
       write(*,*)
     endif
 
+    ! Testing zone
+
+    if(dotest) then
+
+      call dump_test_value('U','UCIS singlet excitation energy',Om_sc(1))
+
+    end if
+
     deallocate(A_sc,Om_sc)
 
   endif
@@ -128,6 +138,14 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
       write(*,*)
     endif
 
+    ! Testing zone
+
+    if(dotest) then
+
+      call dump_test_value('U','UCIS triplet excitation energy',Om_sf(1))
+
+    end if
+
     deallocate(A_sf,Om_sf)
 
   endif
diff --git a/src/GF/GG0F2.f90 b/src/GF/GG0F2.f90
index 4b82ef5..ee03df4 100644
--- a/src/GF/GG0F2.f90
+++ b/src/GF/GG0F2.f90
@@ -1,4 +1,4 @@
-subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
+subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
                  nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform a one-shot second-order Green function calculation
@@ -8,6 +8,8 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: dophBSE
   logical,intent(in)            :: doppBSE
   logical,intent(in)            :: TDA
@@ -39,10 +41,11 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
 
 ! Hello world
 
+
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|     One-shot second-order Green function     |'
-  write(*,*)'************************************************'
+  write(*,*)'********************************'
+  write(*,*)'* Generalized G0F2 Calculation *'
+  write(*,*)'********************************'
   write(*,*)
 
 ! Memory allocation
@@ -80,7 +83,7 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
 
   ! Print results
 
-  call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec)
+  call GMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec)
   call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
 
 ! Perform BSE2 calculation
@@ -113,4 +116,14 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
 
 ! end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('G','GG0F2 correlation energy',Ec)
+    call dump_test_value('G','GG0F2 HOMO energy',eGF(nO))
+    call dump_test_value('G','GG0F2 LUMO energy',eGF(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GF/GGF.f90 b/src/GF/GGF.f90
index 8570bfa..c8064ae 100644
--- a/src/GF/GGF.f90
+++ b/src/GF/GGF.f90
@@ -1,4 +1,4 @@
-subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & 
+subroutine GGF(dotest,doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,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)
 
 ! Green's function module
@@ -8,9 +8,11 @@ subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA
 
 ! Input variables
 
-  logical                       :: doG0F2
-  logical                       :: doevGF2
-  logical                       :: doqsGF2
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doG0F2
+  logical,intent(in)            :: doevGF2
+  logical,intent(in)            :: doqsGF2
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -61,7 +63,7 @@ subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA
   if(doG0F2) then
 
     call wall_time(start_GF)
-    call GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & 
+    call GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & 
                nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GF)
 
@@ -78,7 +80,7 @@ subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA
   if(doevGF2) then
 
     call wall_time(start_GF)
-    call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,linearize,eta,regularize, & 
+    call evGGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,linearize,eta,regularize, & 
                 nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GF)
 
@@ -95,7 +97,7 @@ subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA
   if(doqsGF2) then 
 
     call wall_time(start_GF)
-!   call qsGGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & 
+!   call qsGGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,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_GF)
 
diff --git a/src/GF/G0F2.f90 b/src/GF/RG0F2.f90
similarity index 80%
rename from src/GF/G0F2.f90
rename to src/GF/RG0F2.f90
index daa307a..3c5fd33 100644
--- a/src/GF/G0F2.f90
+++ b/src/GF/RG0F2.f90
@@ -1,5 +1,5 @@
-subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
-                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
+                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform a one-shot second-order Green function calculation
 
@@ -8,6 +8,8 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: dophBSE
   logical,intent(in)            :: doppBSE
   logical,intent(in)            :: TDA
@@ -42,9 +44,9 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|     One-shot second-order Green function     |'
-  write(*,*)'************************************************'
+  write(*,*)'*******************************'
+  write(*,*)'* Restricted G0F2 Calculation *'
+  write(*,*)'*******************************'
   write(*,*)
 
 ! Memory allocation
@@ -82,7 +84,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
 
   ! Print results
 
-  call RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
+  call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
   call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
 
 ! Perform BSE2 calculation
@@ -95,8 +97,8 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2  correlation energy (singlet) =',EcBSE(1)
     write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2  correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2  correlation energy           =',sum(EcBSE(:))
-    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2  total energy                 =',ENuc + ERHF + sum(EcBSE(:))
+    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2  correlation energy           =',sum(EcBSE)
+    write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2  total energy                 =',ENuc + ERHF + sum(EcBSE)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -108,15 +110,27 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
    
     call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
 
+    EcBSE(2) = 3d0*EcBSE(2)
+
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',3d0*EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy           =',EcBSE(1) + 3d0*EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy                 =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',EcBSE(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy           =',sum(EcBSE),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy                 =',ENuc + ERHF + sum(EcBSE),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','RG0F2 correlation energy',Ec)
+    call dump_test_value('R','RG0F2 HOMO energy',eGF(nO))
+    call dump_test_value('R','RG0F2 LUMO energy',eGF(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GF/G0F3.f90 b/src/GF/RG0F3.f90
similarity index 99%
rename from src/GF/G0F3.f90
rename to src/GF/RG0F3.f90
index 90e67a6..4995415 100644
--- a/src/GF/G0F3.f90
+++ b/src/GF/RG0F3.f90
@@ -1,4 +1,4 @@
- subroutine G0F3(renormalization,nBas,nC,nO,nV,nR,V,e0)
+ subroutine RG0F3(dotest,renormalization,nBas,nC,nO,nV,nR,V,e0)
 
 ! Perform third-order Green function calculation in diagonal approximation
 
@@ -7,6 +7,7 @@
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
   integer,intent(in)            :: renormalization
   integer,intent(in)            :: nBas,nC,nO,nV,nR
   double precision,intent(in)   :: e0(nBas),V(nBas,nBas,nBas,nBas)
diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90
index 6ff348e..490dd49 100644
--- a/src/GF/RGF.f90
+++ b/src/GF/RGF.f90
@@ -1,4 +1,4 @@
-subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis,    &
+subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis,    &
                dophBSE,doppBSE,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,11 +10,13 @@ subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
 
 ! Input variables
 
-  logical                       :: doG0F2
-  logical                       :: doevGF2
-  logical                       :: doqsGF2
-  logical                       :: doG0F3
-  logical                       :: doevGF3
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doG0F2
+  logical,intent(in)            :: doevGF2
+  logical,intent(in)            :: doqsGF2
+  logical,intent(in)            :: doG0F3
+  logical,intent(in)            :: doevGF3
 
   integer                       :: renorm
   integer,intent(in)            :: maxSCF
@@ -68,8 +70,8 @@ subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doG0F2) then
 
     call wall_time(start_GF)
-    call G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & 
-              nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & 
+               nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_GF)
 
     t_GF = end_GF - start_GF
@@ -85,9 +87,9 @@ subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doevGF2) then
 
     call wall_time(start_GF)
-    call evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & 
-               singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & 
-               ERI,dipole_int,epsHF)
+    call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & 
+                singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & 
+                ERI,dipole_int,epsHF)
     call wall_time(end_GF)
 
     t_GF = end_GF - start_GF
@@ -103,8 +105,8 @@ subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doqsGF2) then 
 
     call wall_time(start_GF)
-    call qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,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 qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,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_GF)
 
     t_GF = end_GF - start_GF
@@ -120,7 +122,7 @@ subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doG0F3) then
 
     call wall_time(start_GF)
-    call G0F3(renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
+    call RG0F3(dotest,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
     call wall_time(end_GF)
 
     t_GF = end_GF - start_GF
@@ -136,7 +138,7 @@ subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doevGF3) then
 
     call wall_time(start_GF)
-    call evGF3(maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
+    call evRGF3(dotest,maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
     call wall_time(end_GF)
 
     t_GF = end_GF - start_GF
diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90
index 8383c96..a90bfa5 100644
--- a/src/GF/UG0F2.f90
+++ b/src/GF/UG0F2.f90
@@ -1,4 +1,4 @@
-subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
+subroutine UG0F2(dotest,BSE,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,eHF)
 
 ! Perform unrestricted G0W0 calculation
@@ -9,6 +9,8 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regula
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: BSE
   logical,intent(in)            :: TDA
   logical,intent(in)            :: dBSE
@@ -44,18 +46,17 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regula
   double precision,allocatable  :: Z(:,:)
   integer                       :: nS_aa,nS_bb,nS_sc
 
-  double precision,allocatable  :: eGF2lin(:,:)
-  double precision,allocatable  :: eGF2(:,:)
+  double precision,allocatable  :: eGFlin(:,:)
+  double precision,allocatable  :: eGF(:,:)
 
 ! Output variables
 
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|          One-shot G0F2 calculation           |'
-  write(*,*)'|         *** Unrestricted version ***         |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Unrestricted G0F2 Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! TDA 
@@ -71,7 +72,7 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regula
   nS_bb = nS(2)
   nS_sc = nS_aa + nS_bb
 
-  allocate(SigC(nBas,nspin),Z(nBas,nspin),eGF2(nBas,nspin),eGF2lin(nBas,nspin))
+  allocate(SigC(nBas,nspin),Z(nBas,nspin),eGF(nBas,nspin),eGFlin(nBas,nspin))
 
 !---------------------!
 ! Compute self-energy !
@@ -91,14 +92,14 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regula
 ! Solve the quasi-particle equation !
 !-----------------------------------!
 
-  eGF2lin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
+  eGFlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
 
   if(linearize) then 
  
     write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
     write(*,*)
 
-    eGF2(:,:) = eGF2lin(:,:)
+    eGF(:,:) = eGFlin(:,:)
 
   else 
   
@@ -111,11 +112,11 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regula
 
 ! Compute MP2 correlation energy
 
-  call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eGF2,Ec)
+  call UMP2(.false.,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eGF,Ec)
 
 ! Dump results
 
-  call print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
+  call print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF,Ec)
 
 ! Perform BSE calculation
 
@@ -125,4 +126,16 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regula
 
   end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('U','UG0F2 correlation energy',Ec)
+    call dump_test_value('U','UG0F2 HOMOa energy',eGF(nO(1),1))
+    call dump_test_value('U','UG0F2 LUMOa energy',eGF(nO(1)+1,1))
+    call dump_test_value('U','UG0F2 HOMOa energy',eGF(nO(2),2))
+    call dump_test_value('U','UG0F2 LUMOa energy',eGF(nO(2)+1,2))
+
+  end if
+
 end subroutine 
diff --git a/src/GF/UGF.f90 b/src/GF/UGF.f90
index 1aa8357..55bafb9 100644
--- a/src/GF/UGF.f90
+++ b/src/GF/UGF.f90
@@ -1,4 +1,4 @@
-subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis,                      &
+subroutine UGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis,                      &
                dophBSE,doppBSE,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,11 +10,13 @@ subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
 
 ! Input variables
 
-  logical                       :: doG0F2
-  logical                       :: doevGF2
-  logical                       :: doqsGF2
-  logical                       :: doG0F3
-  logical                       :: doevGF3
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: doG0F2
+  logical,intent(in)            :: doevGF2
+  logical,intent(in)            :: doqsGF2
+  logical,intent(in)            :: doG0F3
+  logical,intent(in)            :: doevGF3
 
   integer                       :: renorm
   integer,intent(in)            :: maxSCF
@@ -71,7 +73,7 @@ subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doG0F2) then
 
     call wall_time(start_GF)
-    call UG0F2(dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize, &
+    call UG0F2(dotest,dophBSE,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,epsHF)
     call wall_time(end_GF)
@@ -89,7 +91,7 @@ subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doevGF2) then
 
     call wall_time(start_GF)
-    call evUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,  &
+    call evUGF2(dotest,maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,  &
                 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_GF)
@@ -107,7 +109,7 @@ subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_di
   if(doqsGF2) then 
 
     call wall_time(start_GF)
-    call qsUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, &
+    call qsUGF2(dotest,maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,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)
     call wall_time(end_GF)
diff --git a/src/GF/evGGF2.f90 b/src/GF/evGGF2.f90
index 086c7ff..3fc0c62 100644
--- a/src/GF/evGGF2.f90
+++ b/src/GF/evGGF2.f90
@@ -1,4 +1,4 @@
-subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
+subroutine evGGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
                   linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform eigenvalue self-consistent second-order Green function calculation
@@ -8,6 +8,8 @@ subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: dophBSE
   logical,intent(in)            :: doppBSE
   logical,intent(in)            :: TDA
@@ -50,9 +52,9 @@ subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|  Second-order Green function calculation     |'
-  write(*,*)'************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Generalized evGF2 Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! Memory allocation
@@ -105,7 +107,7 @@ subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
 
     ! Print results
 
-    call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
+    call GMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
     call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
 
     ! DIIS extrapolation
@@ -170,4 +172,14 @@ subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
 
 ! end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('G','evGGF2 correlation energy',Ec)
+    call dump_test_value('G','evGGF2 HOMO energy',eGF(nO))
+    call dump_test_value('G','evGGF2 LUMO energy',eGF(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GF/evGF2.f90 b/src/GF/evRGF2.f90
similarity index 90%
rename from src/GF/evGF2.f90
rename to src/GF/evRGF2.f90
index ebc85da..fd34854 100644
--- a/src/GF/evGF2.f90
+++ b/src/GF/evRGF2.f90
@@ -1,4 +1,4 @@
-subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, &
+subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, &
                  linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform eigenvalue self-consistent second-order Green function calculation
@@ -8,6 +8,8 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: dophBSE
   logical,intent(in)            :: doppBSE
   logical,intent(in)            :: TDA
@@ -51,10 +53,11 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
 
 ! Hello world
 
+
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|  Second-order Green function calculation     |'
-  write(*,*)'************************************************'
+  write(*,*)'********************************'
+  write(*,*)'* Restricted evGF2 Calculation *'
+  write(*,*)'********************************'
   write(*,*)
 
 ! Memory allocation
@@ -107,7 +110,7 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
 
     ! Print results
 
-    call RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
+    call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
     call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
 
     ! DIIS extrapolation
@@ -176,4 +179,14 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
 
   end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','evRGF2 correlation energy',Ec)
+    call dump_test_value('R','evRGF2 HOMO energy',eGF(nO))
+    call dump_test_value('R','evRGF2 LUMO energy',eGF(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GF/evGF3.f90 b/src/GF/evRGF3.f90
similarity index 98%
rename from src/GF/evGF3.f90
rename to src/GF/evRGF3.f90
index 8451ee5..852cbc1 100644
--- a/src/GF/evGF3.f90
+++ b/src/GF/evRGF3.f90
@@ -1,4 +1,4 @@
- subroutine evGF3(maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
+ subroutine evRGF3(dotest,maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
 
 ! Perform third-order Green function calculation in diagonal approximation
 
@@ -7,6 +7,8 @@
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   double precision,intent(in)   :: thresh
   integer,intent(in)            :: maxSCF,max_diis,renormalization
   integer,intent(in)            :: nBas,nC,nO,nV,nR
diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90
index e01066d..d19c5d0 100644
--- a/src/GF/evUGF2.f90
+++ b/src/GF/evUGF2.f90
@@ -1,4 +1,4 @@
-subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
+subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,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,8 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -61,10 +63,11 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
 
 ! Hello world
 
+
   write(*,*)
-  write(*,*)'**************************************************'
-  write(*,*)'| Self-consistent unrestricted evGF2 calculation |'
-  write(*,*)'**************************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Unrestricted G0F2 Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! TDA 
@@ -133,7 +136,7 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
 
     ! Compute MP2 correlation energy
 
-    call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eGF2,Ec)
+    call UMP2(.false.,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eGF2,Ec)
 
     ! Print results
     
diff --git a/src/GF/qsGF2.f90 b/src/GF/qsRGF2.f90
similarity index 93%
rename from src/GF/qsGF2.f90
rename to src/GF/qsRGF2.f90
index 25eb2f1..2f13450 100644
--- a/src/GF/qsGF2.f90
+++ b/src/GF/qsRGF2.f90
@@ -1,4 +1,4 @@
-subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, & 
+subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,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)
 
@@ -9,6 +9,8 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -78,10 +80,11 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
 
 ! Hello world
 
+
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|       Self-consistent qsGF2 calculation      |'
-  write(*,*)'************************************************'
+  write(*,*)'********************************'
+  write(*,*)'* Restricted qsGF2 Calculation *'
+  write(*,*)'********************************'
   write(*,*)
 
 ! Warning 
@@ -216,7 +219,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
 
     ! Correlation energy
 
-    call RMP2(regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec)
+    call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec)
 
     ! Total energy
 
@@ -288,4 +291,14 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
 
   end if
 
+! Testing zone
+
+  if(dotest) then
+
+    call dump_test_value('R','qsRGF2 correlation energy',Ec)
+    call dump_test_value('R','qsRGF2 HOMO energy',eGF(nO))
+    call dump_test_value('R','qsRGF2 LUMO energy',eGF(nO+1))
+
+  end if
+
 end subroutine 
diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90
index 5ba102c..a2c56cb 100644
--- a/src/GF/qsUGF2.f90
+++ b/src/GF/qsUGF2.f90
@@ -1,4 +1,4 @@
-subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, &
+subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, &
                  nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,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,8 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -89,10 +91,11 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
 
 ! Hello world
 
+
   write(*,*)
-  write(*,*)'**************************************************'
-  write(*,*)'| Self-consistent unrestricted qsGF2 calculation |'
-  write(*,*)'**************************************************'
+  write(*,*)'**********************************'
+  write(*,*)'* Unrestricted qsGF2 Calculation *'
+  write(*,*)'**********************************'
   write(*,*)
 
 ! Warning 
@@ -290,7 +293,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
 
     ! Correlation energy
 
-    call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EqsGF2,eGF2,Ec)
+    call UMP2(.false.,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EqsGF2,eGF2,Ec)
 
     ! Total energy
 
diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90
index 7070099..4514402 100644
--- a/src/QuAcK/GQuAcK.f90
+++ b/src/QuAcK/GQuAcK.f90
@@ -233,7 +233,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
   if(doGF) then
 
     call wall_time(start_GF)
-    call GGF(doG0F2,doevGF2,doqsGF2,maxSCF_GF,thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,lin_GF,eta_GF,reg_GF, &
+    call GGF(dotest,doG0F2,doevGF2,doqsGF2,maxSCF_GF,thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,lin_GF,eta_GF,reg_GF, &
              nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
     call wall_time(end_GF)
 
@@ -252,7 +252,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
   if(doGW) then
     
     call wall_time(start_GW)
-    call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS,    & 
+    call GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS,    & 
              dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, & 
              nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
     call wall_time(end_GW)
diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90
index ecdfc5d..f2d1d97 100644
--- a/src/QuAcK/RQuAcK.f90
+++ b/src/QuAcK/RQuAcK.f90
@@ -225,7 +225,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
   if(doCC) then
 
     call wall_time(start_CC)
-    call RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & 
+    call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & 
              maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call wall_time(end_CC)
 
@@ -244,7 +244,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
   if(doCI) then
 
     call wall_time(start_CI)
-    call RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
+    call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
              epsHF,EHF,cHF,S)
     call wall_time(end_CI)
 
@@ -282,7 +282,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
   if(doGF) then
 
     call wall_time(start_GF)
-    call RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & 
+    call RGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & 
              dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF,              &
              nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,            &
              dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
@@ -303,7 +303,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
   if(doGW) then
     
     call wall_time(start_GW)
-    call RGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
+    call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
              exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet,     &
              lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,            &  
              ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
@@ -324,7 +324,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
   if(doGT) then
     
     call wall_time(start_GT)
-    call RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,  &
+    call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,  &
              exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,             &
              lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,                    &
              ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90
index a2c7af4..4cd189f 100644
--- a/src/QuAcK/UQuAcK.f90
+++ b/src/QuAcK/UQuAcK.f90
@@ -259,7 +259,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
   if(doCI) then
 
     call wall_time(start_CI)
-    call UCI(doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, & 
+    call UCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, & 
              ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,EHF,cHF,S)
     call wall_time(end_CI)
 
@@ -297,7 +297,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
   if(doGF) then
 
     call wall_time(start_GF)
-    call UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF,          &
+    call UGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF,          &
              dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF,              &
              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)
@@ -318,7 +318,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
   if(doGW) then
     
     call wall_time(start_GW)
-    call UGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,      &
+    call UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,      &
              exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
              lin_GW,eta_GW,reg_GW,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)
@@ -339,7 +339,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
   if(doGT) then
     
     call wall_time(start_GT)
-    call UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,  &
+    call UGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,  &
              exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,    &
              lin_GT,eta_GT,reg_GT,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)
diff --git a/src/RPA/crRRPA.f90 b/src/RPA/crRRPA.f90
index 6176fad..53bc65c 100644
--- a/src/RPA/crRRPA.f90
+++ b/src/RPA/crRRPA.f90
@@ -105,10 +105,10 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA correlation energy (singlet) =',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA correlation energy (triplet) =',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA correlation energy           =',sum(EcRPA),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA total energy                 =',ENuc + EHF + sum(EcRPA),' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -125,13 +125,19 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcRPA(1)
-    write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy           =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2)
+    write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA correlation energy (singlet) =',EcRPA(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA correlation energy (triplet) =',EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA correlation energy           =',sum(EcRPA),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA total energy                 =',ENuc + EHF + sum(EcRPA),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+  
+    call dump_test_value('R','crRRPA correlation energy',sum(EcRPA))
+  
+  end if
+
 end subroutine 
diff --git a/test/Gtest_ref.dat b/test/Gtest_ref.dat
index 4c0052e..6c5e3d4 100644
--- a/test/Gtest_ref.dat
+++ b/test/Gtest_ref.dat
@@ -14,3 +14,15 @@
   -0.368057033788489
  ppGRPA correlation energy
   -0.092561239023951
+ GG0F2 correlation energy
+  -0.139599095383592
+ GG0F2 HOMO energy
+  -0.398113441401304
+ GG0F2 LUMO energy
+   0.189499108558060
+ evGGF2 correlation energy
+  -0.143423802624327
+ evGGF2 HOMO energy
+  -0.374883361361466
+ evGGF2 LUMO energy
+   0.188029526841004
diff --git a/test/Rtest_ref.dat b/test/Rtest_ref.dat
index 7c20f60..ae36d46 100644
--- a/test/Rtest_ref.dat
+++ b/test/Rtest_ref.dat
@@ -6,13 +6,45 @@
    0.203278954950938
  RHF dipole moment
    0.611349538338893
- ROHF energy
- -85.160473714509976
  RMP2 correlation energy
-  -0.128988144386404
+  -0.128988144318865
+ CCD correlation energy
+  -0.134831644514626
+ DCD correlation energy
+  -0.137864343293459
+ CCSD correlation energy
+  -0.135521483429938
+ drCCD correlation energy
+  -0.138552794834488
+ rCCD correlation energy
+  -0.184028508433731
+ crCCD correlation energy
+   0.116371767570026
+ lCCD correlation energy
+  -0.092561235738883
+ pCCD correlation energy
+  -0.033110922532151
+ RCIS singlet excitation energy
+   0.345631869593838
+ RCIS triplet excitation energy
+   0.310356023571962
  phRRPA correlation energy
-  -0.138552809856833
+  -0.138552809810619
  phRRPAx correlation energy
-  -0.197284981952336
+  -0.197284981858331
+ crRRPA correlation energy
+  -0.161690198064065
  ppRRPA correlation energy
-  -0.092561239071529
+  -0.092561239024299
+ RG0F2 correlation energy
+  -0.139599095383591
+ RG0F2 HOMO energy
+  -0.398113441401301
+ RG0F2 LUMO energy
+   0.189499108558074
+ evRGF2 correlation energy
+  -0.143423802624325
+ evRGF2 HOMO energy
+  -0.374883361361465
+ evRGF2 LUMO energy
+   0.188029526841018
diff --git a/test/Utest_ref.dat b/test/Utest_ref.dat
index 45d2df7..d69b0a9 100644
--- a/test/Utest_ref.dat
+++ b/test/Utest_ref.dat
@@ -12,6 +12,10 @@
    0.611349538338891
  UMP2 correlation energy
   -0.128988144318865
+ UCIS singlet excitation energy
+   0.310356023571988
+ UCIS triplet excitation energy
+   0.310356023571958
  phURPA correlation energy
   -0.138552809810790
  phURPAx correlation energy
diff --git a/test/methods.test b/test/methods.test
index 682b927..5921001 100644
--- a/test/methods.test
+++ b/test/methods.test
@@ -1,20 +1,20 @@
 # RHF UHF GHF ROHF 
-  T   T   T   T   
+  T   F   T   F   
 # MP2 MP3 
   T   T   
 # CCD pCCD DCD CCSD CCSD(T) 
-  F   F    F   F    F
+  T   T    T   T    F
 # drCCD rCCD crCCD lCCD
-  F     F    F     F
+  T     T    T     T
 # CIS CIS(D) CID CISD FCI
-  F   F      F   F    F
+  T   F      F   F    F
 # phRPA phRPAx crRPA ppRPA 
-  T     T      F     T 
+  T     T      T     T 
 # G0F2 evGF2 qsGF2 G0F3 evGF3
-  F    F     F     F    F
+  T    T     F     F    F
 # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
   F    F    F    F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
   F      F      F      F      F      F
 # Rtest Utest Gtest
-  T     T     T
+  T     F     T
diff --git a/test/options.test b/test/options.test
index 979a4b1..980f0d4 100644
--- a/test/options.test
+++ b/test/options.test
@@ -4,8 +4,8 @@
       F
 # CC: maxSCF thresh     DIIS
       64     0.0000001  5
-# spin: TDA spin_conserved spin_flip 
-        F   T               T
+# spin: TDA singlet triplet
+        F   T       T
 # 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