diff --git a/input/methods b/input/methods
index 73b5c29..682b927 100644
--- a/input/methods
+++ b/input/methods
@@ -1,7 +1,7 @@
 # RHF UHF GHF ROHF 
   T   T   T   T   
 # MP2 MP3 
-  T   F   
+  T   T   
 # CCD pCCD DCD CCSD CCSD(T) 
   F   F    F   F    F
 # drCCD rCCD crCCD lCCD
@@ -9,7 +9,7 @@
 # CIS CIS(D) CID CISD FCI
   F   F      F   F    F
 # phRPA phRPAx crRPA ppRPA 
-  F     F      F     F 
+  T     T      F     T 
 # G0F2 evGF2 qsGF2 G0F3 evGF3
   F    F     F     F    F
 # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
@@ -17,4 +17,4 @@
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
   F      F      F      F      F      F
 # Rtest Utest Gtest
-  F     F     F
+  T     T     T
diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90
index e2131bf..de30055 100644
--- a/src/HF/GHF.f90
+++ b/src/HF/GHF.f90
@@ -1,4 +1,4 @@
-subroutine GHF(doGtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
+subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
                nBas,nBas2,nO,Ov,T,V,Hc,ERI,dipole_int,Or,EHF,e,c,P)
 
 ! Perform unrestricted Hartree-Fock calculation
@@ -8,7 +8,7 @@ subroutine GHF(doGtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
 
 ! Input variables
 
-  logical,intent(in)            :: doGtest
+  logical,intent(in)            :: dotest
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -298,7 +298,7 @@ subroutine GHF(doGtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
 
 ! Print test values
 
-  if(doGtest) then
+  if(dotest) then
 
     call dump_test_value('G','GHF energy',EHF)
 
diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90
index 477f823..1867cdf 100644
--- a/src/HF/RHF.f90
+++ b/src/HF/RHF.f90
@@ -1,4 +1,4 @@
-subroutine RHF(doRtest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
+subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
                nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P)
 
 ! Perform restricted Hartree-Fock calculation
@@ -8,7 +8,7 @@ subroutine RHF(doRtest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 
 ! Input variables
 
-  logical,intent(in)            :: doRtest
+  logical,intent(in)            :: dotest
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -203,7 +203,7 @@ subroutine RHF(doRtest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 
 ! Print test values
 
-  if(doRtest) then
+  if(dotest) then
  
     call dump_test_value('R','RHF energy',EHF)
 
diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90
index 2606c25..a84ef20 100644
--- a/src/HF/ROHF.f90
+++ b/src/HF/ROHF.f90
@@ -1,4 +1,4 @@
-subroutine ROHF(doRtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
+subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
                 nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,Ptot)
 
 ! Perform restricted open-shell Hartree-Fock calculation
@@ -8,7 +8,7 @@ subroutine ROHF(doRtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,Z
 
 ! Input variables
 
-  logical,intent(in)            :: doRtest
+  logical,intent(in)            :: dotest
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -245,7 +245,7 @@ subroutine ROHF(doRtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,Z
 
 ! Print test values
 
-  if(doRtest) then
+  if(dotest) then
 
     call dump_test_value('R','ROHF energy',EHF)
 
diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90
index 0d6791c..bc89cf1 100644
--- a/src/HF/UHF.f90
+++ b/src/HF/UHF.f90
@@ -1,4 +1,4 @@
-subroutine UHF(doUtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
+subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
                nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P)
 
 ! Perform unrestricted Hartree-Fock calculation
@@ -8,7 +8,7 @@ subroutine UHF(doUtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
 
 ! Input variables
 
-  logical,intent(in)            :: doUtest
+  logical,intent(in)            :: dotest
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -256,7 +256,7 @@ subroutine UHF(doUtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
 
 ! Print test values
 
-  if(doUtest) then
+  if(dotest) then
 
     call dump_test_value('U','UHF energy',EHF)
 
diff --git a/src/MP/GMP.f90 b/src/MP/GMP.f90
index eeb239e..7249803 100644
--- a/src/MP/GMP.f90
+++ b/src/MP/GMP.f90
@@ -1,4 +1,4 @@
-subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+subroutine GMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
 
 ! Moller-Plesset module
 
@@ -7,6 +7,8 @@ subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: doMP2
   logical,intent(in)            :: doMP3
 
@@ -35,7 +37,7 @@ subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
   if(doMP2) then    
        
     call wall_time(start_MP)
-    call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec)
+    call GMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
diff --git a/src/MP/GMP2.f90 b/src/MP/GMP2.f90
index 48be3b7..9558629 100644
--- a/src/MP/GMP2.f90
+++ b/src/MP/GMP2.f90
@@ -1,4 +1,4 @@
-subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
+subroutine GMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF,EcMP2)
 
 ! Perform second-order Moller-Plesset calculation with and without regularizers
 
@@ -6,6 +6,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: regularize
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -13,8 +15,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: EHF
-  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: EGHF
+  double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
 ! Local variables
@@ -38,9 +40,9 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|    Moller-Plesset second-order calculation   |'
-  write(*,*)'************************************************'
+  write(*,*)'*******************************'
+  write(*,*)'* Generalized MP2 Calculation |'
+  write(*,*)'*******************************'
   write(*,*)
 
 !---------------------------------------------!
@@ -70,7 +72,7 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
       do a=nO+1,nBas-nR
         do b=nO+1,nBas-nR
 
-          Dijab = e(a) + e(b) - e(i) - e(j) 
+          Dijab = eHF(a) + eHF(b) - eHF(i) - eHF(j) 
 
           ! Second-order ring diagram
 
@@ -114,8 +116,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
   write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2d
   write(*,'(A32,1X,F16.10)') ' Exchange part           = ',E2x
   write(*,'(A32)')           '---------------------------'
-  write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + EcMP2
-  write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + EcMP2
+  write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EGHF + EcMP2
+  write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EGHF + EcMP2
   write(*,'(A32)')           '---------------------------'
   write(*,*)
 
@@ -133,8 +135,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
     write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2ds
     write(*,'(A32,1X,F16.10)') ' Exchange part           = ',E2xs
     write(*,'(A32)')           '---------------------------'
-    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + EcsMP2
-    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + EcsMP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EGHF + EcsMP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EGHF + EcsMP2
     write(*,'(A32)')           '---------------------------'
     write(*,*)
 
@@ -150,8 +152,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
     write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2ds2
     write(*,'(A32,1X,F16.10)') ' Exchange part           = ',E2xs2
     write(*,'(A32)')           '---------------------------'
-    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + Ecs2MP2
-    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + Ecs2MP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EGHF + Ecs2MP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EGHF + Ecs2MP2
     write(*,'(A32)')           '---------------------------'
     write(*,*)
 
@@ -167,11 +169,18 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
     write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2dk
     write(*,'(A32,1X,F16.10)') ' Exchange part           = ',-E2xk
     write(*,'(A32)')           '---------------------------'
-    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + EckMP2
-    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + EckMP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EGHF + EckMP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EGHF + EckMP2
     write(*,'(A32)')           '---------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('G','GMP2 correlation energy',EcMP2)
+
+  end if
+
+
 end subroutine 
diff --git a/src/MP/GMP3.f90 b/src/MP/GMP3.f90
index 5f2e494..70ced8b 100644
--- a/src/MP/GMP3.f90
+++ b/src/MP/GMP3.f90
@@ -38,9 +38,9 @@ subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|    Moller-Plesset third-order calculation    |'
-  write(*,*)'************************************************'
+  write(*,*)'*******************************'
+  write(*,*)'* Generalized MP3 Calculation *'
+  write(*,*)'*******************************'
   write(*,*)
 
 ! Antysymmetrize ERIs
diff --git a/src/MP/RMP.f90 b/src/MP/RMP.f90
index d3c8a06..0a88427 100644
--- a/src/MP/RMP.f90
+++ b/src/MP/RMP.f90
@@ -1,4 +1,4 @@
-subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
+subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! Moller-Plesset module
 
@@ -7,6 +7,8 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: doMP2
   logical,intent(in)            :: doMP3
 
@@ -35,7 +37,7 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
   if(doMP2) then    
        
     call wall_time(start_MP)
-    call RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,Ec)
+    call RMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,Ec)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
diff --git a/src/MP/RMP2.f90 b/src/MP/RMP2.f90
index c6a10d5..f944599 100644
--- a/src/MP/RMP2.f90
+++ b/src/MP/RMP2.f90
@@ -1,4 +1,4 @@
-subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
+subroutine RMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
 
 ! Perform second-order Moller-Plesset calculation with and without regularizers
 
@@ -6,6 +6,8 @@ subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: regularize
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -37,9 +39,9 @@ subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|    Moller-Plesset second-order calculation   |'
-  write(*,*)'************************************************'
+  write(*,*)'******************************'
+  write(*,*)'* Restricted MP2 Calculation *'
+  write(*,*)'******************************'
   write(*,*)
 
 !---------------------------------------------!
@@ -171,4 +173,10 @@ subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('R','RMP2 correlation energy',EcMP2)
+
+  end if
+
 end subroutine 
diff --git a/src/MP/RMP3.f90 b/src/MP/RMP3.f90
index c5a8b72..fc7b817 100644
--- a/src/MP/RMP3.f90
+++ b/src/MP/RMP3.f90
@@ -45,9 +45,9 @@ subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e)
 ! Hello world
 
   write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|    Moller-Plesset third-order calculation    |'
-  write(*,*)'************************************************'
+  write(*,*)'******************************'
+  write(*,*)'* Restricted MP3 Calculation *'
+  write(*,*)'******************************'
   write(*,*)
 
 ! Spatial to spin orbitals
diff --git a/src/MP/UMP.f90 b/src/MP/UMP.f90
index 6d1dc46..21f77f7 100644
--- a/src/MP/UMP.f90
+++ b/src/MP/UMP.f90
@@ -1,4 +1,4 @@
-subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
+subroutine UMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
 
 ! Moller-Plesset module
 
@@ -7,6 +7,8 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: doMP2
   logical,intent(in)            :: doMP3
 
@@ -37,7 +39,7 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb
   if(doMP2) then    
        
     call wall_time(start_MP)
-    call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF,Ec)
+    call UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF,Ec)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
diff --git a/src/MP/UMP2.f90 b/src/MP/UMP2.f90
index d7a361c..3c19874 100644
--- a/src/MP/UMP2.f90
+++ b/src/MP/UMP2.f90
@@ -1,24 +1,25 @@
-subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
+subroutine UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EUHF,eHF,Ec)
 
 ! Perform unrestricted second-order Moller-Plesset calculation
 
   implicit none
   include 'parameters.h'
 
-
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)
   integer,intent(in)            :: nO(nspin)
   integer,intent(in)            :: nV(nspin)
   integer,intent(in)            :: nR(nspin)
   double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: EUHF
   double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: eHF(nBas,nspin)
 
 ! Local variables
 
@@ -37,9 +38,9 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
 ! Hello world
 
   write(*,*)
-  write(*,*)'********************************************************'
-  write(*,*)'| Unrestricted second-order Moller-Plesset calculation |'
-  write(*,*)'********************************************************'
+  write(*,*)'********************************'
+  write(*,*)'* Unrestricted MP2 Calculation *'
+  write(*,*)'********************************'
   write(*,*)
 
 !---------------------!
@@ -60,7 +61,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
       do j=nC(ket)+1,nO(ket)
         do b=nO(ket)+1,nBas-nR(ket)
 
-          eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket) 
+          eps = eHF(i,bra) + eHF(j,ket) - eHF(a,bra) - eHF(b,ket) 
          
           Edaa = Edaa + 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,a,b)/eps
           Exaa = Exaa - 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,b,a)/eps
@@ -88,7 +89,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
       do j=nC(ket)+1,nO(ket)
         do b=nO(ket)+1,nBas-nR(ket)
 
-          eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket) 
+          eps = eHF(i,bra) + eHF(j,ket) - eHF(a,bra) - eHF(b,ket) 
          
           Edab = Edab + ERI_ab(i,j,a,b)*ERI_ab(i,j,a,b)/eps
 
@@ -114,7 +115,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
       do j=nC(ket)+1,nO(ket)
         do b=nO(ket)+1,nBas-nR(ket)
 
-          eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket) 
+          eps = eHF(i,bra) + eHF(j,ket) - eHF(a,bra) - eHF(b,ket) 
          
           Edbb = Edbb + 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,a,b)/eps
           Exbb = Exbb - 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,b,a)/eps
@@ -152,9 +153,15 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
   write(*,'(A32,1X,F16.10)') '   alpha-beta           = ',Exab
   write(*,'(A32,1X,F16.10)') '    beta-beta           = ',Exbb
   write(*,'(A32)')           '--------------------------'
-  write(*,'(A32,1X,F16.10)') ' MP2 electronic  energy = ',       EHF + sum(Ec(:))
-  write(*,'(A32,1X,F16.10)') ' MP2 total       energy = ',ENuc + EHF + sum(Ec(:))
+  write(*,'(A32,1X,F16.10)') ' MP2 electronic  energy = ',       EUHF + sum(Ec(:))
+  write(*,'(A32,1X,F16.10)') ' MP2 total       energy = ',ENuc + EUHF + sum(Ec(:))
   write(*,'(A32)')           '--------------------------'
   write(*,*)
 
+  if(dotest) then
+
+    call dump_test_value('U','UMP2 correlation energy',sum(Ec))
+
+  end if
+
 end subroutine 
diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90
index 35ba25e..7070099 100644
--- a/src/QuAcK/GQuAcK.f90
+++ b/src/QuAcK/GQuAcK.f90
@@ -1,4 +1,4 @@
-subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,   &
+subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,   &
                   doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,                          &
                   nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
                   maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP,    &
@@ -9,7 +9,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
   implicit none
   include 'parameters.h'
 
-  logical,intent(in)            :: doGtest
+  logical,intent(in)            :: dotest
 
   logical,intent(in)            :: doGHF
   logical,intent(in)            :: dostab
@@ -103,7 +103,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
   if(doGHF) then
 
     call wall_time(start_HF)
-    call GHF(doGtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
+    call GHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
@@ -196,7 +196,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
   if(doMP) then
 
     call wall_time(start_MP)
-    call GMP(doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
+    call GMP(dotest,doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
@@ -214,7 +214,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
   if(doRPA) then
 
     call wall_time(start_RPA)
-    call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, & 
+    call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, & 
               ERI_MO,dipole_int_MO,epsHF,cHF,S)
     call wall_time(end_RPA)
 
diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90
index d6993a6..1975b07 100644
--- a/src/QuAcK/RQuAcK.f90
+++ b/src/QuAcK/RQuAcK.f90
@@ -1,4 +1,4 @@
-subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,              &   
+subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,              &   
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & 
                   doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,            &
                   doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,      & 
@@ -13,7 +13,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
   implicit none
   include 'parameters.h'
 
-  logical,intent(in)            :: doRtest
+  logical,intent(in)            :: dotest
 
   logical,intent(in)            :: doRHF,doROHF
   logical,intent(in)            :: dostab
@@ -117,11 +117,12 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
   if(doRHF) then
 
     call wall_time(start_HF)
-    call RHF(doRtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
+    call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
     t_HF = end_HF - start_HF
+    write(*,*)
     write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
     write(*,*)
 
@@ -130,7 +131,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
   if(doROHF) then
 
     call wall_time(start_HF)
-    call ROHF(doRtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
+    call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
               nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
@@ -206,7 +207,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
   if(doMP) then
 
     call wall_time(start_MP)
-    call RMP(doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
+    call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
@@ -263,7 +264,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
   if(doRPA) then
 
     call wall_time(start_RPA)
-    call RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet,    & 
+    call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet,    & 
               nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
     call wall_time(end_RPA)
 
diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90
index 4927ec8..a2c7af4 100644
--- a/src/QuAcK/UQuAcK.f90
+++ b/src/QuAcK/UQuAcK.f90
@@ -1,4 +1,4 @@
-subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                     &
+subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                     &
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
                   doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,            &
                   doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,      &
@@ -11,7 +11,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
   implicit none
   include 'parameters.h'
 
-  logical,intent(in)            :: doUtest
+  logical,intent(in)            :: dotest
 
   logical,intent(in)            :: doUHF
   logical,intent(in)            :: dostab
@@ -116,7 +116,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
   if(doUHF) then
 
     call wall_time(start_HF)
-    call UHF(doUtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
+    call UHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
@@ -223,7 +223,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
   if(doMP) then
 
     call wall_time(start_MP)
-    call UMP(doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
+    call UMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
@@ -278,7 +278,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
   if(doRPA) then
 
     call wall_time(start_RPA)
-    call URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,    & 
+    call URPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,    & 
               nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
     call wall_time(end_RPA)
 
diff --git a/src/RPA/GRPA.f90 b/src/RPA/GRPA.f90
index ce90a38..745fa05 100644
--- a/src/RPA/GRPA.f90
+++ b/src/RPA/GRPA.f90
@@ -1,4 +1,4 @@
-subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
 
 ! Random-phase approximation module
 
@@ -7,9 +7,11 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
 
 ! Input variables
 
-  logical                       :: dophRPA
-  logical                       :: dophRPAx
-  logical                       :: doppRPA
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: dophRPA
+  logical,intent(in)            :: dophRPAx
+  logical,intent(in)            :: doppRPA
 
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
@@ -37,7 +39,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
   if(dophRPA) then
 
     call wall_time(start_RPA)
-    call phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -53,7 +55,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
   if(dophRPAx) then
 
     call wall_time(start_RPA)
-    call phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
 
     call wall_time(end_RPA)
 
@@ -70,7 +72,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
   if(doppRPA) then
 
     call wall_time(start_RPA)
-    call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
+    call ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
diff --git a/src/RPA/RRPA.f90 b/src/RPA/RRPA.f90
index 4c36d77..2fc1706 100644
--- a/src/RPA/RRPA.f90
+++ b/src/RPA/RRPA.f90
@@ -1,4 +1,4 @@
-subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &  
+subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &  
                 nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF,cHF,S)
 
 ! Random-phase approximation module
@@ -8,10 +8,12 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
 
 ! Input variables
 
-  logical                       :: dophRPA
-  logical                       :: dophRPAx
-  logical                       :: docrRPA
-  logical                       :: doppRPA
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: dophRPA
+  logical,intent(in)            :: dophRPAx
+  logical,intent(in)            :: docrRPA
+  logical,intent(in)            :: doppRPA
 
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
@@ -43,7 +45,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
   if(dophRPA) then
 
     call wall_time(start_RPA)
-    call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -59,7 +61,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
   if(dophRPAx) then
 
     call wall_time(start_RPA)
-    call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -75,7 +77,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
   if(docrRPA) then
 
     call wall_time(start_RPA)
-    call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -91,7 +93,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
   if(doppRPA) then
 
     call wall_time(start_RPA)
-    call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
+    call ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
diff --git a/src/RPA/URPA.f90 b/src/RPA/URPA.f90
index a312009..a8b3e96 100644
--- a/src/RPA/URPA.f90
+++ b/src/RPA/URPA.f90
@@ -1,4 +1,4 @@
-subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &  
+subroutine URPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &  
                 nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
 
 ! Random-phase approximation module
@@ -8,10 +8,12 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
 
 ! Input variables
 
-  logical                       :: dophRPA
-  logical                       :: dophRPAx
-  logical                       :: docrRPA
-  logical                       :: doppRPA
+  logical,intent(in)            :: dotest
+
+  logical,intent(in)            :: dophRPA
+  logical,intent(in)            :: dophRPAx
+  logical,intent(in)            :: docrRPA
+  logical,intent(in)            :: doppRPA
 
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
@@ -46,7 +48,7 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
   if(dophRPA) then
 
     call wall_time(start_RPA)
-    call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
+    call phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
                 ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
     call wall_time(end_RPA)
 
@@ -63,7 +65,7 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
   if(dophRPAx) then
 
     call wall_time(start_RPA)
-    call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
+    call phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
                  ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
     call wall_time(end_RPA)
 
@@ -96,7 +98,7 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
   if(doppRPA) then
 
     call wall_time(start_RPA)
-    call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
+    call ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRRPA.f90
similarity index 96%
rename from src/RPA/crRPA.f90
rename to src/RPA/crRRPA.f90
index b5283d7..6176fad 100644
--- a/src/RPA/crRPA.f90
+++ b/src/RPA/crRRPA.f90
@@ -1,4 +1,4 @@
-subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Crossed-ring channel of the random phase approximation 
 
@@ -8,6 +8,8 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
diff --git a/src/RPA/phGRPA.f90 b/src/RPA/phGRPA.f90
index 3ae3e6e..486abea 100644
--- a/src/RPA/phGRPA.f90
+++ b/src/RPA/phGRPA.f90
@@ -1,4 +1,4 @@
-subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Perform a direct random phase approximation calculation
 
@@ -8,6 +8,8 @@ subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -36,9 +38,9 @@ subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 ! Hello world
 
   write(*,*)
-  write(*,*)'***********************************************'
-  write(*,*)'|  Random-phase approximation calculation     |'
-  write(*,*)'***********************************************'
+  write(*,*)'**********************************'
+  write(*,*)'* Generalized ph-RPA Calculation |'
+  write(*,*)'**********************************'
   write(*,*)
 
 ! TDA 
@@ -70,9 +72,15 @@ subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy           =',EcRPA
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy                 =',ENuc + EHF + EcRPA
+  write(*,'(2X,A50,F20.10)') 'Tr@phGRPA correlation energy           = ',EcRPA
+  write(*,'(2X,A50,F20.10)') 'Tr@phGRPA total energy                 = ',ENuc + EHF + EcRPA
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
+  if(dotest) then
+
+    call dump_test_value('G','phGRPA corrlation energy',EcRPA)
+
+  end if
+
 end subroutine 
diff --git a/src/RPA/phGRPAx.f90 b/src/RPA/phGRPAx.f90
index 7fc1a7a..d1c38a3 100644
--- a/src/RPA/phGRPAx.f90
+++ b/src/RPA/phGRPAx.f90
@@ -1,4 +1,4 @@
-subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF)
 
@@ -8,6 +8,8 @@ subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -36,9 +38,9 @@ subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 ! Hello world
 
   write(*,*)
-  write(*,*)'***********************************************************'
-  write(*,*)'|  Random phase approximation calculation with exchange   |'
-  write(*,*)'***********************************************************'
+  write(*,*)'************************************'
+  write(*,*)'* Generalized ph-RPAx Calculation  *'
+  write(*,*)'************************************'
   write(*,*)
 
 ! TDA 
@@ -71,9 +73,15 @@ subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy           =',EcRPA
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy                 =',ENuc + EHF + EcRPA
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPAx correlation energy           = ',EcRPA,' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPAx total energy                 = ',ENuc + EHF + EcRPA,' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
+  if(dotest) then
+
+    call dump_test_value('G','phGRPAx correlation energy',EcRPA)
+
+  end if
+
 end subroutine 
diff --git a/src/RPA/phRPA.f90 b/src/RPA/phRRPA.f90
similarity index 77%
rename from src/RPA/phRPA.f90
rename to src/RPA/phRRPA.f90
index 14cb769..8d15258 100644
--- a/src/RPA/phRPA.f90
+++ b/src/RPA/phRRPA.f90
@@ -1,4 +1,4 @@
-subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Perform a direct random phase approximation calculation
 
@@ -8,6 +8,8 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
@@ -102,10 +104,10 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA correlation energy (singlet) = ',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA correlation energy (triplet) = ',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA correlation energy           = ',sum(EcRPA),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA total energy                 = ',ENuc + EHF + sum(EcRPA),' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -124,13 +126,19 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (singlet) =',EcRPA(1)
-    write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (triplet) =',EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy           =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2)
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA correlation energy (singlet) = ',EcRPA(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA correlation energy (triplet) = ',EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA correlation energy           = ',sum(EcRPA),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA total energy                 = ',ENuc + EHF + sum(EcRPA),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('R','phRRPA correlation energy',sum(EcRPA))
+
+  end if
+
 end subroutine 
diff --git a/src/RPA/phRPAx.f90 b/src/RPA/phRRPAx.f90
similarity index 77%
rename from src/RPA/phRPAx.f90
rename to src/RPA/phRRPAx.f90
index 1ae61b7..3423ddf 100644
--- a/src/RPA/phRPAx.f90
+++ b/src/RPA/phRRPAx.f90
@@ -1,4 +1,4 @@
-subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF)
 
@@ -8,6 +8,8 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
@@ -103,10 +105,10 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx correlation energy (singlet) = ',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx correlation energy (triplet) = ',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx correlation energy           = ',sum(EcRPA),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx total energy                 = ',ENuc + EHF + sum(EcRPA),' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -127,13 +129,19 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (singlet) =',EcRPA(1)
-    write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (triplet) =',EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy           =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@phRPAx total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2)
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx correlation energy (singlet) = ',EcRPA(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx correlation energy (triplet) = ',EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx correlation energy           = ',sum(EcRPA),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx total energy                 = ',ENuc + EHF + sum(EcRPA),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('R','phRRPAx correlation energy',sum(EcRPA))
+
+  end if
+
 end subroutine 
diff --git a/src/RPA/phURPA.f90 b/src/RPA/phURPA.f90
index 356f12f..c373929 100644
--- a/src/RPA/phURPA.f90
+++ b/src/RPA/phURPA.f90
@@ -1,5 +1,5 @@
-subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & 
-                  ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S)
+subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & 
+                  ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,c,S)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
 
@@ -9,6 +9,8 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
@@ -22,7 +24,7 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
   integer,intent(in)            :: nS(nspin)
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: EUHF
-  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: eHF(nBas,nspin)
   double precision,intent(in)   :: c(nBas,nBas,nspin)
   double precision,intent(in)   :: S(nBas,nBas)
   double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
@@ -81,7 +83,7 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
 
     allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
 
-    call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
+    call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
     if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
 
     call phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
@@ -104,9 +106,9 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
     nSb = (nO(2) - nC(2))*(nV(1) - nR(1))
     nSt = nSa + nSb
 
-    allocate(Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
+    allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
 
-    call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
+    call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
     if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
 
     call phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
@@ -126,10 +128,10 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-conserved) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-flip)      =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy                  =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@URPA total energy                        =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA correlation energy (spin-conserved) = ',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA correlation energy (spin-flip)      = ',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA correlation energy                  = ',sum(EcRPA),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA total energy                        = ',ENuc + EUHF + sum(EcRPA),' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -143,7 +145,7 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
     write(*,*)
 
     call phUACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,spin_conserved,spin_flip, &
-                  nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcRPA)
+                  nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,EcRPA)
 
     if(exchange_kernel) then
    
@@ -154,13 +156,19 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-conserved) =',EcRPA(1)
-    write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-flip)      =',EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy                  =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@URPA total energy                        =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA correlation energy (spin-conserved) = ',EcRPA(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA correlation energy (spin-flip)      = ',EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA correlation energy                  = ',sum(EcRPA),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA total energy                        = ',ENuc + EUHF + sum(EcRPA),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('U','phURPA correlation energy',sum(EcRPA))
+
+  end if
+
 end subroutine 
diff --git a/src/RPA/phURPAx.f90 b/src/RPA/phURPAx.f90
index 0e4c568..1814a35 100644
--- a/src/RPA/phURPAx.f90
+++ b/src/RPA/phURPAx.f90
@@ -1,4 +1,4 @@
-subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & 
+subroutine phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & 
                    ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
@@ -9,6 +9,8 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
@@ -104,7 +106,7 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
     nSb = (nO(2) - nC(2))*(nV(1) - nR(1))
     nSt = nSa + nSb
 
-    allocate(Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
+    allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
 
     call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
     if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
@@ -137,10 +139,10 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-conserved) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-flip)      =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy                  =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@URPAx total energy                        =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx correlation energy (spin-conserved) = ',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx correlation energy (spin-flip)      = ',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx correlation energy                  = ',sum(EcRPA),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx total energy                        = ',ENuc + EUHF + sum(EcRPA),' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -158,13 +160,19 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-conserved) =',EcRPA(1)
-    write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-flip)      =',EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy                  =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A50,F20.10)') 'AC@URPAx total energy                        =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
+    write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx correlation energy (spin-conserved) =',EcRPA(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx correlation energy (spin-flip)      =',EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx correlation energy                  =',sum(EcRPA),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx total energy                        =',ENuc + EUHF + sum(EcRPA),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('U','phURPAx correlation energy',sum(EcRPA))
+
+  end if
+
 end subroutine 
diff --git a/src/RPA/ppGRPA.f90 b/src/RPA/ppGRPA.f90
index 45c2ad8..1136f99 100644
--- a/src/RPA/ppGRPA.f90
+++ b/src/RPA/ppGRPA.f90
@@ -1,4 +1,4 @@
-subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
+subroutine ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 
 ! Perform ppGRPA calculation
 
@@ -7,6 +7,8 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   integer,intent(in)            :: nBas
@@ -40,9 +42,9 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 ! Hello world
 
   write(*,*)
-  write(*,*)'****************************************'
-  write(*,*)'|  particle-particle GRPA calculation  |'
-  write(*,*)'****************************************'
+  write(*,*)'**********************************'
+  write(*,*)'* Generalized pp-RPA Calculation *'
+  write(*,*)'**********************************'
   write(*,*)
 
 ! Initialization
@@ -70,8 +72,8 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy           =',EcRPA
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy                 =',ENuc + EHF + EcRPA
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppGRPA correlation energy           = ',EcRPA,' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppGRPA total energy                 = ',ENuc + EHF + EcRPA,' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -88,13 +90,17 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 
 !   write(*,*)
 !   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcRPA(1),' au'
-!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcRPA(2),' au'
-!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy           =',EcRPA(1) + EcRPA(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy           =',EcRPA,' au'
 !   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
 !   write(*,*)'-------------------------------------------------------------------------------'
 !   write(*,*)
 
 ! end if
 
+  if(dotest) then
+  
+    call dump_test_value('G','ppGRPA correlation energy',EcRPA)
+
+  end if
+
 end subroutine
diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRRPA.f90
similarity index 77%
rename from src/RPA/ppRPA.f90
rename to src/RPA/ppRRPA.f90
index 9da9357..660076c 100644
--- a/src/RPA/ppRPA.f90
+++ b/src/RPA/ppRRPA.f90
@@ -1,4 +1,4 @@
-subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
+subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 
 ! Perform ppRPA calculation
 
@@ -7,6 +7,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: singlet
@@ -42,9 +44,9 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
 ! Hello world
 
   write(*,*)
-  write(*,*)'****************************************'
-  write(*,*)'|  particle-particle RPA calculation   |'
-  write(*,*)'****************************************'
+  write(*,*)'*********************************'
+  write(*,*)'* Restricted pp-RPA Calculation *'
+  write(*,*)'*********************************'
   write(*,*)
 
 ! Initialization
@@ -115,12 +117,14 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
 
   endif
 
+  EcRPA(2) = 3d0*EcRPA(2)
+
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy           =',EcRPA(1) + 3d0*EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy                 =',ENuc + EHF + EcRPA(1) + 3d0*EcRPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy (singlet) = ',EcRPA(1),'au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy (triplet) = ',EcRPA(2),'au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy           = ',sum(EcRPA),'au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA total energy                 = ',ENuc + EHF + sum(EcRPA),'au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -137,13 +141,19 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcRPA(1),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcRPA(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy           =',EcRPA(1) + EcRPA(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy                 =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy (singlet) = ',EcRPA(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy (triplet) = ',EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy           = ',EcRPA(1) + EcRPA(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA total energy                 = ',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
+  if(dotest) then
+
+    call dump_test_value('R','ppRRPA correlation energy',sum(EcRPA))
+
+  end if
+
 end subroutine
diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90
index 55aa283..6380a65 100644
--- a/src/RPA/ppURPA.f90
+++ b/src/RPA/ppURPA.f90
@@ -1,4 +1,4 @@
-subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,e)
+subroutine ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,e)
 
 ! Perform unrestricted pp-RPA calculations
 
@@ -7,6 +7,8 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
 
 ! Input variables
 
+  logical,intent(in)            :: dotest
+
   logical,intent(in)            :: TDA
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: spin_conserved
@@ -35,20 +37,20 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
   double precision,allocatable  :: X2sc(:,:),X2sf(:,:)
   double precision,allocatable  :: Y2sc(:,:),Y2sf(:,:)
 
-  double precision              :: Ec_ppURPA(nspin)
+  double precision              :: EcRPA(nspin)
   double precision              :: EcAC(nspin)
 
 ! Hello world
 
   write(*,*)
-  write(*,*)'****************************************'
-  write(*,*)'|  particle-particle URPA calculation  |'
-  write(*,*)'****************************************'
+  write(*,*)'***********************************'
+  write(*,*)'* Unrestricted pp-RPA Calculation *'
+  write(*,*)'***********************************'
   write(*,*)
 
 ! Initialization
 
-  Ec_ppURPA(:) = 0d0
+  EcRPA(:) = 0d0
   EcAC(:)   = 0d0
 
 !alpha-beta block  
@@ -70,7 +72,7 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
   call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
              nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1sc,X1sc,Y1sc, &
-             Om2sc,X2sc,Y2sc,Ec_ppURPA(ispin))
+             Om2sc,X2sc,Y2sc,EcRPA(ispin))
     
   call print_excitation_energies('ppRPA@UHF (N+2)',iblock,nP_sc,Om1sc)
   call print_excitation_energies('ppRPA@UHF (N-2)',iblock,nH_sc,Om2sc)
@@ -96,7 +98,7 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
   call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
              nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf, &
-             Om2sf,X2sf,Y2sf,Ec_ppURPA(ispin))
+             Om2sf,X2sf,Y2sf,EcRPA(ispin))
 
   call print_excitation_energies('ppRPA@UHF (N+2)',iblock,nP_sf,Om1sf)
   call print_excitation_energies('ppRPA@UHF (N-2)',iblock,nH_sf,Om2sf)
@@ -116,17 +118,19 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
   call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,&
              nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,&
              ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf,&
-             Om2sf,X2sf,Y2sf,Ec_ppURPA(ispin))
+             Om2sf,X2sf,Y2sf,EcRPA(ispin))
 
   call print_excitation_energies('ppRPA@UHF (N+2)',iblock,nP_sf,Om1sf)
   call print_excitation_energies('ppRPA@UHF (N-2)',iblock,nH_sf,Om2sf)
 
+  EcRPA(2) = 3d0*EcRPA(2)
+
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-conserved) =',Ec_ppURPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-flip)      =',3d0*Ec_ppURPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy                  =',Ec_ppURPA(1) + 3d0*Ec_ppURPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy                        =',ENuc + EUHF + Ec_ppURPA(1) + 3d0*Ec_ppURPA(2)
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA correlation energy (spin-conserved) = ',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA correlation energy (spin-flip)      = ',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA correlation energy                  = ',sum(EcRPA),' au'
+  write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA total energy                        = ',ENuc + EUHF + sum(EcRPA),' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -152,4 +156,10 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
 
 ! end if
 
+  if(dotest) then
+
+    call dump_test_value('U','ppURPA correlation energy',sum(EcRPA))
+
+  end if
+
 end subroutine 
diff --git a/src/test/run_test.f90 b/src/test/run_test.f90
index 6bd6b0a..669cbf0 100755
--- a/src/test/run_test.f90
+++ b/src/test/run_test.f90
@@ -10,6 +10,8 @@ subroutine run_test(doRtest,doUtest,doGtest)
  
 ! Local variables
 
+  double precision              :: start_test     ,end_test       ,t_test
+
 ! Output variables
 
   if(doRtest) then 
@@ -19,7 +21,12 @@ subroutine run_test(doRtest,doUtest,doGtest)
     write(*,*) '****************************************'
     write(*,*) 
 
+    call wall_time(start_test)
     call check_test_value('R')
+    call wall_time(end_test)
+
+    t_test = end_test - start_test
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for test of restricted branch = ',t_test,' seconds'
 
     write(*,*) 
     write(*,*) '**************************'
@@ -36,7 +43,12 @@ subroutine run_test(doRtest,doUtest,doGtest)
     write(*,*) '******************************************'
     write(*,*) 
 
+    call wall_time(start_test)
     call check_test_value('U')
+    call wall_time(end_test)
+
+    t_test = end_test - start_test
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for test of unrestricted branch = ',t_test,' seconds'
 
     write(*,*) 
     write(*,*) '****************************'
@@ -53,7 +65,12 @@ subroutine run_test(doRtest,doUtest,doGtest)
     write(*,*) '*****************************************'
     write(*,*) 
 
+    call wall_time(start_test)
     call check_test_value('G')
+    call wall_time(end_test)
+
+    t_test = end_test - start_test
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for test of generalized branch = ',t_test,' seconds'
 
     write(*,*) 
     write(*,*) '***************************'
diff --git a/test/Gtest_ref.dat b/test/Gtest_ref.dat
index 94106f6..f711ac2 100644
--- a/test/Gtest_ref.dat
+++ b/test/Gtest_ref.dat
@@ -1,2 +1,10 @@
- # GHF energy
+ GHF energy
  -85.160473883160876
+ GMP2 correlation energy
+  -0.128988144318866
+ phGRPA corrlation energy
+  -0.138552809810392
+ phGRPAx correlation energy
+  -0.368057033788489
+ ppGRPA correlation energy
+  -0.092561239023951
diff --git a/test/Rtest_ref.dat b/test/Rtest_ref.dat
index 370de7d..9e0aee7 100644
--- a/test/Rtest_ref.dat
+++ b/test/Rtest_ref.dat
@@ -2,3 +2,11 @@
  -85.160473883160876
  ROHF energy
  -85.160473714509976
+ RMP2 correlation energy
+  -0.128988144386404
+ phRRPA correlation energy
+  -0.138552809856833
+ phRRPAx correlation energy
+  -0.197284981952336
+ ppRRPA correlation energy
+  -0.092561239071529
diff --git a/test/Utest_ref.dat b/test/Utest_ref.dat
index 38f3b77..9ab06b9 100644
--- a/test/Utest_ref.dat
+++ b/test/Utest_ref.dat
@@ -1,2 +1,10 @@
- # UHF energy
+ UHF energy
  -85.160473883160819
+ UMP2 correlation energy
+  -0.128988144318865
+ phURPA correlation energy
+  -0.138552809810790
+ phURPAx correlation energy
+  -0.197284981858218
+ ppURPA correlation energy
+  -0.103998858975444
diff --git a/test/methods.test b/test/methods.test
index 9fbdc3c..682b927 100644
--- a/test/methods.test
+++ b/test/methods.test
@@ -1,7 +1,7 @@
 # RHF UHF GHF ROHF 
   T   T   T   T   
 # MP2 MP3 
-  F   F   
+  T   T   
 # CCD pCCD DCD CCSD CCSD(T) 
   F   F    F   F    F
 # drCCD rCCD crCCD lCCD
@@ -9,7 +9,7 @@
 # CIS CIS(D) CID CISD FCI
   F   F      F   F    F
 # phRPA phRPAx crRPA ppRPA 
-  F     F      F     F 
+  T     T      F     T 
 # G0F2 evGF2 qsGF2 G0F3 evGF3
   F    F     F     F    F
 # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW