diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90
index de30055..d8e24d7 100644
--- a/src/HF/GHF.f90
+++ b/src/HF/GHF.f90
@@ -1,5 +1,5 @@
 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)
+               nBas,nBas2,nO,Ov,T,V,Hc,ERI,dipole_int,Or,EGHF,eHF,c,P)
 
 ! Perform unrestricted Hartree-Fock calculation
 
@@ -66,8 +66,8 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
 ! Output variables
 
-  double precision,intent(out)  :: EHF
-  double precision,intent(out)  :: e(nBas2)
+  double precision,intent(out)  :: EGHF
+  double precision,intent(out)  :: eHF(nBas2)
   double precision,intent(inout):: C(nBas2,nBas2)
   double precision,intent(out)  :: P(nBas2,nBas2)
 
@@ -202,7 +202,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 !  Diagonalize Fock matrix to get eigenvectors and eigenvalues
 
     Cp(:,:) = Fp(:,:)
-    call diagonalize_matrix(nBas2,Cp,e)
+    call diagonalize_matrix(nBas2,Cp,eHF)
     
 !   Back-transform eigenvectors in non-orthogonal basis
 
@@ -261,12 +261,12 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
 !   Total energy
 
-    EHF = ET + EV + EJ + EK
+    EGHF = ET + EV + EJ + EK
 
 !   Dump results
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & 
-      '|',nSCF,'|',EHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
+      '|',nSCF,'|',EGHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
  
   end do
   write(*,*)'-----------------------------------------------------------------------------'
@@ -290,17 +290,20 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
 ! Compute dipole moments
 
-  call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
+  call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
 
 ! Compute final GHF energy
 
-  call print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
+  call print_GHF(nBas,nBas2,nO,eHF,C,P,ENuc,ET,EV,EJ,EK,EGHF,dipole)
 
 ! Print test values
 
   if(dotest) then
 
-    call dump_test_value('G','GHF energy',EHF)
+    call dump_test_value('G','GHF energy',EGHF)
+    call dump_test_value('G','GHF HOMO energy',eHF(nO))
+    call dump_test_value('G','GHF LUMO energy',eHF(nO+1))
+    call dump_test_value('G','GHF dipole moment',norm2(dipole))
 
   end if
 
diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90
index 1867cdf..f07e4bd 100644
--- a/src/HF/RHF.f90
+++ b/src/HF/RHF.f90
@@ -1,5 +1,5 @@
 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)
+               nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,eHF,c,P)
 
 ! Perform restricted Hartree-Fock calculation
 
@@ -56,8 +56,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 
 ! Output variables
 
-  double precision,intent(out)  :: EHF
-  double precision,intent(out)  :: e(nBas)
+  double precision,intent(out)  :: ERHF
+  double precision,intent(out)  :: eHF(nBas)
   double precision,intent(inout):: c(nBas,nBas)
   double precision,intent(out)  :: P(nBas,nBas)
 
@@ -138,7 +138,7 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 
     Fp = matmul(transpose(X),matmul(F,X))
     cp(:,:) = Fp(:,:)
-    call diagonalize_matrix(nBas,cp,e)
+    call diagonalize_matrix(nBas,cp,eHF)
     c = matmul(X,cp)
 
 !   Density matrix
@@ -147,15 +147,15 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
   
 !   Compute HF energy
 
-    EHF = trace_matrix(nBas,matmul(P,Hc))      &
-        + 0.5d0*trace_matrix(nBas,matmul(P,J)) &
-        + 0.25d0*trace_matrix(nBas,matmul(P,K))
+    ERHF = trace_matrix(nBas,matmul(P,Hc))      &
+         + 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+         + 0.25d0*trace_matrix(nBas,matmul(P,K))
 
 !   Compute HOMO-LUMO gap
 
     if(nBas > nO) then 
 
-      Gap = e(nO+1) - e(nO)
+      Gap = eHF(nO+1) - eHF(nO)
 
     else
 
@@ -166,7 +166,7 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 !  Dump results
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & 
-      '|',nSCF,'|',EHF+ENuc,'|',Conv,'|',Gap,'|'
+      '|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|'
 
   end do
   write(*,*)'----------------------------------------------------'
@@ -194,18 +194,21 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
   EV = trace_matrix(nBas,matmul(P,V))
   EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
   EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
-  EHF = ET + EV + EJ + EK
+  ERHF = ET + EV + EJ + EK
 
 ! Compute dipole moments
 
   call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
-  call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,EHF,dipole)
+  call print_RHF(nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
 
 ! Print test values
 
   if(dotest) then
  
-    call dump_test_value('R','RHF energy',EHF)
+    call dump_test_value('R','RHF energy',ERHF)
+    call dump_test_value('R','RHF HOMO energy',eHF(nO))
+    call dump_test_value('R','RHF LUMO energy',eHF(nO+1))
+    call dump_test_value('R','RHF dipole moment',norm2(dipole))
 
   end if
 
diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90
index bc89cf1..1e30d3f 100644
--- a/src/HF/UHF.f90
+++ b/src/HF/UHF.f90
@@ -1,5 +1,5 @@
 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)
+               nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,eHF,c,P)
 
 ! Perform unrestricted Hartree-Fock calculation
 
@@ -59,8 +59,8 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
 ! Output variables
 
-  double precision,intent(out)  :: EHF
-  double precision,intent(out)  :: e(nBas,nspin)
+  double precision,intent(out)  :: EUHF
+  double precision,intent(out)  :: eHF(nBas,nspin)
   double precision,intent(inout):: c(nBas,nBas,nspin)
   double precision,intent(out)  :: P(nBas,nBas,nspin)
 
@@ -172,7 +172,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
     cp(:,:,:) = Fp(:,:,:)
     do ispin=1,nspin
-      call diagonalize_matrix(nBas,cp(:,:,ispin),e(:,ispin))
+      call diagonalize_matrix(nBas,cp(:,:,ispin),eHF(:,ispin))
     end do
     
 !   Back-transform eigenvectors in non-orthogonal basis
@@ -221,12 +221,12 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
 !   Total energy
 
-    EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
+    EUHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
 
 !   Dump results
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & 
-      '|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|'
+      '|',nSCF,'|',EUHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|'
  
   end do
   write(*,*)'----------------------------------------------------------'
@@ -251,14 +251,19 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 ! Compute final UHF energy
 
   call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
-  call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
+  call print_UHF(nBas,nO,S,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
 
 
 ! Print test values
 
   if(dotest) then
 
-    call dump_test_value('U','UHF energy',EHF)
+    call dump_test_value('U','UHF energy',EUHF)
+    call dump_test_value('U','UHF HOMOa energy',eHF(nO(1),1))
+    call dump_test_value('U','UHF HOMOb energy',eHF(nO(2),2))
+    call dump_test_value('U','UHF LUMOa energy',eHF(nO(1)+1,1))
+    call dump_test_value('U','UHF LUMOb energy',eHF(nO(2)+1,2))
+    call dump_test_value('U','UHF dipole moment',norm2(dipole))
 
   end if
 
diff --git a/src/test/check_test_value.f90 b/src/test/check_test_value.f90
index f906e34..2ee9bde 100755
--- a/src/test/check_test_value.f90
+++ b/src/test/check_test_value.f90
@@ -51,9 +51,9 @@ subroutine check_test_value(branch)
     read(12,'(F20.15)',end=12) reference
 
     if(abs(value-reference) < cutoff) then
-      answer = '..... [SUCCESS]'
+      answer = '.......... :-)'
     else
-      answer = '..... [FAILED] '
+      answer = '.......... :-( '
       failed = .true.
     end if
     write(*,'(1X,A1,1X,A30,1X,A1,1X,3F15.10,1X,A1,1X,A15,1X,A1)') & 
diff --git a/test/Gtest_ref.dat b/test/Gtest_ref.dat
index f711ac2..4c0052e 100644
--- a/test/Gtest_ref.dat
+++ b/test/Gtest_ref.dat
@@ -1,5 +1,11 @@
  GHF energy
  -85.160473883160876
+ GHF HOMO energy
+  -0.501365804897696
+ GHF LUMO energy
+   0.203278954950924
+ GHF dipole moment
+   0.714967673390535
  GMP2 correlation energy
   -0.128988144318866
  phGRPA corrlation energy
diff --git a/test/Rtest_ref.dat b/test/Rtest_ref.dat
index 9e0aee7..7c20f60 100644
--- a/test/Rtest_ref.dat
+++ b/test/Rtest_ref.dat
@@ -1,5 +1,11 @@
  RHF energy
  -85.160473883160876
+ RHF HOMO energy
+  -0.501365804897693
+ RHF LUMO energy
+   0.203278954950938
+ RHF dipole moment
+   0.611349538338893
  ROHF energy
  -85.160473714509976
  RMP2 correlation energy
diff --git a/test/Utest_ref.dat b/test/Utest_ref.dat
index 9ab06b9..45d2df7 100644
--- a/test/Utest_ref.dat
+++ b/test/Utest_ref.dat
@@ -1,5 +1,15 @@
  UHF energy
  -85.160473883160819
+ UHF HOMOa energy
+  -0.501365804897699
+ UHF HOMOb energy
+  -0.501365804897702
+ UHF LUMOa energy
+   0.203278954950949
+ UHF LUMOb energy
+   0.203278954950948
+ UHF dipole moment
+   0.611349538338891
  UMP2 correlation energy
   -0.128988144318865
  phURPA correlation energy