diff --git a/src/CC/DEA_EOM_CCD_2p.f90 b/src/CC/DEA_EOM_CCD_2p.f90
index 36c0b04..335947a 100644
--- a/src/CC/DEA_EOM_CCD_2p.f90
+++ b/src/CC/DEA_EOM_CCD_2p.f90
@@ -111,6 +111,6 @@ subroutine DEA_EOM_CCD_2p(nC,nO,nV,nR,eV,OOVV,VVVV,t)
 
 ! Dump results
 
-  call print_excitation_energies('DEA-EOM-CCD',3,nVV,Om)
+  call print_excitation_energies('DEA-EOM-CCD','spinorbital',nVV,Om)
 
 end subroutine 
diff --git a/src/CC/DIP_EOM_CCD_2h.f90 b/src/CC/DIP_EOM_CCD_2h.f90
index 28472a4..2e9290a 100644
--- a/src/CC/DIP_EOM_CCD_2h.f90
+++ b/src/CC/DIP_EOM_CCD_2h.f90
@@ -111,6 +111,6 @@ subroutine DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t)
 
 ! Dump results
 
-  call print_excitation_energies('DIP-EOM-CCD',3,nOO,Om)
+  call print_excitation_energies('DIP-EOM-CCD','spinorbital',nOO,Om)
 
 end subroutine 
diff --git a/src/CC/EE_EOM_CCD_1h1p.f90 b/src/CC/EE_EOM_CCD_1h1p.f90
index 55ce7ec..d50ed49 100644
--- a/src/CC/EE_EOM_CCD_1h1p.f90
+++ b/src/CC/EE_EOM_CCD_1h1p.f90
@@ -143,7 +143,7 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
     call quick_sort(Om,order,nS)
     call set_order(Z,order,nS,nS)
 
-    call print_excitation_energies('EE-EOM-CCD',3,nS,Om)
+    call print_excitation_energies('EE-EOM-CCD','spinorbital',nS,Om)
 
   end if
 
diff --git a/src/CI/RCIS.f90 b/src/CI/RCIS.f90
index 54f0add..5bdaeca 100644
--- a/src/CI/RCIS.f90
+++ b/src/CI/RCIS.f90
@@ -57,7 +57,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
     endif
 
     call diagonalize_matrix(nS,A,Om)
-    call print_excitation_energies('CIS@RHF',ispin,nS,Om)
+    call print_excitation_energies('CIS@RHF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,transpose(A),transpose(A))
  
     if(dump_trans) then
@@ -93,7 +93,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
     endif
  
     call diagonalize_matrix(nS,A,Om)
-    call print_excitation_energies('CIS@RHF',ispin,nS,Om)
+    call print_excitation_energies('CIS@RHF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,transpose(A),transpose(A))
 
     if(dump_trans) then
diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90
index 148f9e0..407d92b 100644
--- a/src/CI/UCIS.f90
+++ b/src/CI/UCIS.f90
@@ -80,7 +80,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI
 
     call diagonalize_matrix(nS_sc,A_sc,Om_sc)
     A_sc(:,:) = transpose(A_sc)
-    call print_excitation_energies('CIS@UHF',5,nS_sc,Om_sc)
+    call print_excitation_energies('CIS@UHF','spin-conserved',nS_sc,Om_sc)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
                                   cHF,S,Om_sc,A_sc,A_sc)
  
@@ -128,7 +128,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI
 
     call diagonalize_matrix(nS_sf,A_sf,Om_sf)
     A_sf(:,:) = transpose(A_sf)
-    call print_excitation_energies('CIS@UHF',6,nS_sf,Om_sf)
+    call print_excitation_energies('CIS@UHF','spin-flip',nS_sf,Om_sf)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
                                   cHF,S,Om_sf,A_sf,A_sf)
  
diff --git a/src/GF/GF2_phBSE2.f90 b/src/GF/GF2_phBSE2.f90
index 97ec667..a3e8ad8 100644
--- a/src/GF/GF2_phBSE2.f90
+++ b/src/GF/GF2_phBSE2.f90
@@ -68,7 +68,7 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
     ! Compute phBSE2@GF2 excitation energies
 
     call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
-    call print_excitation_energies('phBSE2@GF2',ispin,nS,OmBSE)
+    call print_excitation_energies('phBSE2@GF2','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
     ! Compute dynamic correction for BSE via perturbation theory
@@ -101,7 +101,7 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
     ! Compute phBSE2@GF2 excitation energies
 
     call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
-    call print_excitation_energies('phBSE2@GF2',ispin,nS,OmBSE)
+    call print_excitation_energies('phBSE2@GF2','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
     ! Compute dynamic correction for BSE via perturbation theory
diff --git a/src/GF/GGF2_phBSE2.f90 b/src/GF/GGF2_phBSE2.f90
index 471d723..d101ce2 100644
--- a/src/GF/GGF2_phBSE2.f90
+++ b/src/GF/GGF2_phBSE2.f90
@@ -60,7 +60,7 @@ subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,
   ! Compute phBSE2@GF2 excitation energies
 
   call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
-  call print_excitation_energies('phBSE2@GGF2',ispin,nS,OmBSE)
+  call print_excitation_energies('phBSE2@GGF2','spinorbital',nS,OmBSE)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
   ! Compute dynamic correction for BSE via perturbation theory
diff --git a/src/GT/GTpp_phBSE.f90 b/src/GT/GTpp_phBSE.f90
index bf72191..a0eac9d 100644
--- a/src/GT/GTpp_phBSE.f90
+++ b/src/GT/GTpp_phBSE.f90
@@ -139,7 +139,7 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
     call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
-    call print_excitation_energies('phBSE@GTpp',ispin,nS,OmBSE)
+    call print_excitation_energies('phBSE@GTpp','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     ! Compute dynamic correction for BSE via renormalized perturbation theory 
@@ -169,7 +169,7 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
     call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
-    call print_excitation_energies('phBSE@GTpp',ispin,nS,OmBSE)
+    call print_excitation_energies('phBSE@GTpp','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     ! Compute dynamic correction for BSE via renormalized perturbation theory 
diff --git a/src/GT/RG0T0eh.f90 b/src/GT/RG0T0eh.f90
index c721acf..1b96ad7 100644
--- a/src/GT/RG0T0eh.f90
+++ b/src/GT/RG0T0eh.f90
@@ -113,7 +113,7 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
 
   call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
-  if(print_T) call print_excitation_energies('phRPA@HF',ispin,nS,Om)
+  if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
 
 !--------------------------!
 ! Compute spectral weights !
diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90
index 179e3c8..2c4ce09 100644
--- a/src/GT/RG0T0pp.f90
+++ b/src/GT/RG0T0pp.f90
@@ -128,8 +128,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   deallocate(Bpp,Cpp,Dpp)
 
-  call print_excitation_energies('ppRPA@HF (N+2)',iblock,nVVs,Om1s(:))
-  call print_excitation_energies('ppRPA@HF (N-2)',iblock,nOOs,Om2s(:))
+  call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s(:))
+  call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s(:))
 
 !----------------------------------------------
 ! alpha-alpha block
@@ -151,8 +151,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   deallocate(Bpp,Cpp,Dpp)
 
-  call print_excitation_energies('ppRPA (N+2)',iblock,nVVt,Om1t)
-  call print_excitation_energies('ppRPA (N-2)',iblock,nOOt,Om2t)
+  call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)
+  call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOt,Om2t)
 
 !----------------------------------------------
 ! Compute excitation densities
diff --git a/src/GT/UG0T0pp.f90 b/src/GT/UG0T0pp.f90
index a66e4de..c499980 100644
--- a/src/GT/UG0T0pp.f90
+++ b/src/GT/UG0T0pp.f90
@@ -116,8 +116,8 @@ subroutine UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,
   call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) 
 
-  call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPab,Om1ab(:))
-  call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHab,Om2ab(:))
+  call print_excitation_energies('ppRPA@UHF','2p (alpha-beta)',nPab,Om1ab(:))
+  call print_excitation_energies('ppRPA@UHF','2h (alpha-beta)',nHab,Om2ab(:))
  
 !----------------------------------------------
 ! alpha-alpha block
@@ -131,8 +131,8 @@ subroutine UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,
   call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))  
   
-  call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPaa,Om1aa(:))
-  call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHaa,Om2aa(:))
+  call print_excitation_energies('ppRPA@UHF','2p (alpha-alpha)',nPaa,Om1aa(:))
+  call print_excitation_energies('ppRPA@UHF','2h (alpha-alpha)',nHaa,Om2aa(:))
 
 !----------------------------------------------
 ! beta-beta block
@@ -146,8 +146,8 @@ subroutine UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,
   call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, &
              ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
 
-  call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPbb,Om1bb(:))
-  call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHbb,Om2bb(:))
+  call print_excitation_energies('ppRPA@UHF','2p (beta-beta)',nPbb,Om1bb(:))
+  call print_excitation_energies('ppRPA@UHF','2h (beta-beta)',nHbb,Om2bb(:))
 
 !----------------------------------------------
 ! Compute T-matrix version of the self-energy 
diff --git a/src/GT/evRGTeh.f90 b/src/GT/evRGTeh.f90
index 415c100..1e6cb92 100644
--- a/src/GT/evRGTeh.f90
+++ b/src/GT/evRGTeh.f90
@@ -127,7 +127,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
 
     call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
-    if(print_T) call print_excitation_energies('phRPA@GTeh',ispin,nS,Om)
+    if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
 
    ! Compute spectral weights
 
diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90
index 37d9736..5b57f45 100644
--- a/src/GT/qsRGTeh.f90
+++ b/src/GT/qsRGTeh.f90
@@ -178,7 +178,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
 
     call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
-    if(print_T) call print_excitation_energies('phRPA@qsGTeh',ispin,nS,Om)
+    if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
 
     ! Compute correlation part of the self-energy 
 
diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90
index ccadf34..972dce8 100644
--- a/src/GW/GG0W0.f90
+++ b/src/GW/GG0W0.f90
@@ -103,7 +103,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
-  if(print_W) call print_excitation_energies('phRPA@GHF',ispin,nS,Om)
+  if(print_W) call print_excitation_energies('phRPA@GHF','spinorbital',nS,Om)
 
 !--------------------------!
 ! Compute spectral weights !
diff --git a/src/GW/GGW_phBSE.f90 b/src/GW/GGW_phBSE.f90
index d11c168..b1b9479 100644
--- a/src/GW/GGW_phBSE.f90
+++ b/src/GW/GGW_phBSE.f90
@@ -109,7 +109,7 @@ subroutine GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,di
 
   call phLR(TDA,nS,Aph,Bph,EcBSE,OmBSE,XpY_BSE,XmY_BSE)
 
-  call print_excitation_energies('phBSE@GGW',ispin,nS,OmBSE)
+  call print_excitation_energies('phBSE@GGW','spinorbital',nS,OmBSE)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     !----------------------------------------------------!
diff --git a/src/GW/GW_phBSE.f90 b/src/GW/GW_phBSE.f90
index 152ad4d..0cf4e18 100644
--- a/src/GW/GW_phBSE.f90
+++ b/src/GW/GW_phBSE.f90
@@ -117,7 +117,7 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
 
     call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
-    call print_excitation_energies('phBSE@GW',ispin,nS,OmBSE)
+    call print_excitation_energies('phBSE@RGW','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     !----------------------------------------------------!
@@ -149,7 +149,7 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
 
     call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
-    call print_excitation_energies('phBSE@GW',ispin,nS,OmBSE)
+    call print_excitation_energies('phBSE@RGW','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     !-------------------------------------------------
diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90
index 73c9dd2..19f2cb9 100644
--- a/src/GW/RG0W0.f90
+++ b/src/GW/RG0W0.f90
@@ -106,7 +106,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
-  if(print_W) call print_excitation_energies('phRPA@RHF',ispin,nS,Om)
+  if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
 
 !--------------------------!
 ! Compute spectral weights !
diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90
index 64e4454..e0b6167 100644
--- a/src/GW/SRG_qsGW.f90
+++ b/src/GW/SRG_qsGW.f90
@@ -195,7 +195,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
 
     tlr = tlr + tlr2 -tlr1
 
-    if(print_W) call print_excitation_energies('phRPA@SRG-qsGW',ispin,nS,Om)
+    if(print_W) call print_excitation_energies('phRPA@RGW','singlet',nS,Om)
 
     ! Compute correlation part of the self-energy 
 
diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90
index 7d2ae15..b6b1f60 100644
--- a/src/GW/UG0W0.f90
+++ b/src/GW/UG0W0.f90
@@ -115,7 +115,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
   call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
   
-  if(print_W) call print_excitation_energies('phRPA@UHF',5,nSt,Om)
+  if(print_W) call print_excitation_energies('phRPA@UHF','spin-conserved',nSt,Om)
 
 !----------------------!
 ! Excitation densities !
diff --git a/src/GW/UGW_phBSE.f90 b/src/GW/UGW_phBSE.f90
index 0770beb..cf97a2b 100644
--- a/src/GW/UGW_phBSE.f90
+++ b/src/GW/UGW_phBSE.f90
@@ -100,7 +100,7 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,  &
 
     call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
                eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
-    call print_excitation_energies('phBSE@UGW',5,nS_sc,OmBSE_sc)
+    call print_excitation_energies('phBSE@UGW','spin-conserved',nS_sc,OmBSE_sc)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
                                   cW,S,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
 
@@ -136,7 +136,7 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,  &
                eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),       & 
                OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
 
-    call print_excitation_energies('phBSE@UGW',6,nS_sf,OmBSE_sf)
+    call print_excitation_energies('phBSE@UGW','spin-flip',nS_sf,OmBSE_sf)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
                                   cW,S,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
 
diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90
index cfe3fa7..5678676 100644
--- a/src/GW/qsGGW.f90
+++ b/src/GW/qsGGW.f90
@@ -257,7 +257,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
 
     call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-    if(print_W) call print_excitation_energies('phRPA@qsGGW',ispin,nS,Om)
+    if(print_W) call print_excitation_energies('phRPA@GGW','spinorbital',nS,Om)
 
 !   Compute correlation part of the self-energy 
 
diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90
index 3ca58b6..9fc39f8 100644
--- a/src/GW/qsRGW.f90
+++ b/src/GW/qsRGW.f90
@@ -180,7 +180,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
 
     call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-    if(print_W) call print_excitation_energies('phRPA@qsRGW',ispin,nS,Om)
+    if(print_W) call print_excitation_energies('phRPA@RGW','singlet',nS,Om)
 
     ! Compute correlation part of the self-energy 
 
diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90
index 4845bdc..28a4140 100644
--- a/src/HF/print_GHF.f90
+++ b/src/HF/print_GHF.f90
@@ -130,8 +130,8 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
   write(*,'(A33,1X,F16.6)')     ' <Sz>                = ',Sz
   write(*,'(A50)')           '---------------------------------------'
   write(*,'(A33,1X,F16.6)')     ' <Sx^2+Sy^2>         = ',S2 - Sz2
-  write(*,'(A33,1X,F16.6)')     ' <Sz**2>             = ',Sz2
-  write(*,'(A33,1X,F16.6)')     ' <S**2>              = ',S2
+  write(*,'(A33,1X,F16.6)')     ' <Sz^2>              = ',Sz2
+  write(*,'(A33,1X,F16.6)')     ' <S^2>               = ',S2
   write(*,'(A50)')           '---------------------------------------'
   write(*,'(A36)')           ' Dipole moment (Debye)    '
   write(*,'(10X,4A10)')      'X','Y','Z','Tot.'
diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90
index 3d7d5da..3a06d31 100644
--- a/src/HF/print_RHF.f90
+++ b/src/HF/print_RHF.f90
@@ -60,8 +60,8 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
   write(*,'(A33,1X,F16.6,A3)')  ' HF LUMO      energy = ',eHF(LUMO)*HaToeV,' eV'
   write(*,'(A33,1X,F16.6,A3)')  ' HF HOMO-LUMO gap    = ',Gap*HaToeV,' eV'
   write(*,'(A50)')           '---------------------------------------'
-  write(*,'(A33,1X,F16.6)')     '  S                 = ',S
-  write(*,'(A33,1X,F16.6)')     ' <S**2>             = ',S2
+  write(*,'(A33,1X,F16.6)')     ' <Sz>                = ',S
+  write(*,'(A33,1X,F16.6)')     ' <S^2>               = ',S2
   write(*,'(A50)')           '---------------------------------------'
   write(*,'(A36)')           ' Dipole moment (Debye)    '
   write(*,'(10X,4A10)')      'X','Y','Z','Tot.'
diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90
index 140bfa0..9d29535 100644
--- a/src/HF/print_ROHF.f90
+++ b/src/HF/print_ROHF.f90
@@ -53,42 +53,42 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
   write(*,'(A60)')              '-------------------------------------------------'
   write(*,'(A40)')              ' Summary              '
   write(*,'(A60)')              '-------------------------------------------------'
-  write(*,'(A40,1X,F16.10,A3)') ' One-electron    energy: ',sum(ET(:))  + sum(EV(:)),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' One-electron a  energy: ',ET(1) + EV(1),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' One-electron b  energy: ',ET(2) + EV(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Kinetic         energy: ',sum(ET(:)),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Kinetic      a  energy: ',ET(1),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Kinetic      b  energy: ',ET(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Potential       energy: ',sum(EV(:)),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Potential    a  energy: ',EV(1),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Potential    b  energy: ',EV(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron    energy = ',sum(ET)  + sum(EV),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron a  energy = ',ET(1) + EV(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron b  energy = ',ET(2) + EV(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic         energy = ',sum(ET),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic      a  energy = ',ET(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic      b  energy = ',ET(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential       energy = ',sum(EV),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential    a  energy = ',EV(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential    b  energy = ',EV(2),' au'
   write(*,'(A60)')              '-------------------------------------------------'
-  write(*,'(A40,1X,F16.10,A3)') ' Two-electron    energy: ',sum(EJ(:)) + sum(Ex(:)),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Hartree         energy: ',sum(EJ(:)),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Hartree      aa energy: ',EJ(1),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Hartree      ab energy: ',EJ(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Hartree      bb energy: ',EJ(3),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy: ',sum(Ex(:)),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Exchange     a  energy: ',Ex(1),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Exchange     b  energy: ',Ex(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron    energy = ',sum(EJ) + sum(Ex),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy = ',EJ(1) + Ex(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy = ',EJ(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy = ',EJ(3) + Ex(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Hartree         energy = ',sum(EJ),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Hartree      aa energy = ',EJ(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Hartree      ab energy = ',EJ(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Hartree      bb energy = ',EJ(3),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy = ',sum(Ex),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange     a  energy = ',Ex(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange     b  energy = ',Ex(2),' au'
   write(*,'(A60)')              '-------------------------------------------------'
-  write(*,'(A40,1X,F16.10,A3)') ' Electronic      energy: ',EROHF,' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Nuclear      repulsion: ',ENuc,' au'
-  write(*,'(A40,1X,F16.10,A3)') ' ROHF            energy: ',EROHF + ENuc,' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Electronic      energy = ',EROHF,' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Nuclear      repulsion = ',ENuc,' au'
+  write(*,'(A40,1X,F16.10,A3)') ' ROHF            energy = ',EROHF + ENuc,' au'
   write(*,'(A60)')              '-------------------------------------------------'
-  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMO a   energy:',HOMO(1)*HatoeV,' eV'
-  write(*,'(A40,1X,F16.6,A3)')  ' ROHF LUMO a   energy:',LUMO(1)*HatoeV,' eV'
-  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV'
+  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMO a     energy = ',HOMO(1)*HatoeV,' eV'
+  write(*,'(A40,1X,F16.6,A3)')  ' ROHF LUMO a     energy = ',LUMO(1)*HatoeV,' eV'
+  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMOa-LUMOa   gap = ',Gap(1)*HatoeV,' eV'
   write(*,'(A60)')              '-------------------------------------------------'
-  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMO b   energy:',HOMO(2)*HatoeV,' eV'
-  write(*,'(A40,1X,F16.6,A3)')  ' ROHF LUMO b   energy:',LUMO(2)*HatoeV,' eV'
-  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMOb-LUMOb gap:',Gap(2)*HatoeV,' eV'
+  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMO b     energy = ',HOMO(2)*HatoeV,' eV'
+  write(*,'(A40,1X,F16.6,A3)')  ' ROHF LUMO b     energy = ',LUMO(2)*HatoeV,' eV'
+  write(*,'(A40,1X,F16.6,A3)')  ' ROHF HOMOb-LUMOb   gap = ',Gap(2)*HatoeV,' eV'
   write(*,'(A60)')              '-------------------------------------------------'
-  write(*,'(A40,1X,F16.6)')     '  S                  :',2d0*S       + 1d0
-  write(*,'(A40,1X,F16.6)')     ' <S**2>              :',S2
+  write(*,'(A40,1X,F16.6)')     ' <Sz>                   = ',S
+  write(*,'(A40,1X,F16.6)')     ' <S^2>                  = ',S2
   write(*,'(A60)')              '-------------------------------------------------'
   write(*,'(A45)')              ' Dipole moment (Debye)    '
   write(*,'(19X,4A10)')         'X','Y','Z','Tot.'
diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90
index 1cd3ff7..466667f 100644
--- a/src/HF/print_UHF.f90
+++ b/src/HF/print_UHF.f90
@@ -61,25 +61,25 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
   write(*,'(A60)')              '---------------------------------------------'
   write(*,'(A40)')              ' Summary                  '
   write(*,'(A60)')              '---------------------------------------------'
-  write(*,'(A40,1X,F16.10,A3)') ' One-electron    energy = ',sum(ET(:))  + sum(EV(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron    energy = ',sum(ET)  + sum(EV),' au'
   write(*,'(A40,1X,F16.10,A3)') ' One-electron a  energy = ',ET(1) + EV(1),' au'
   write(*,'(A40,1X,F16.10,A3)') ' One-electron b  energy = ',ET(2) + EV(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Kinetic         energy = ',sum(ET(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic         energy = ',sum(ET),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Kinetic      a  energy = ',ET(1),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Kinetic      b  energy = ',ET(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Potential       energy = ',sum(EV(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential       energy = ',sum(EV),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Potential    a  energy = ',EV(1),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Potential    b  energy = ',EV(2),' au'
   write(*,'(A60)')              '---------------------------------------------'
-  write(*,'(A40,1X,F16.10,A3)') ' Two-electron    energy = ',sum(EJ(:)) + sum(Ex(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron    energy = ',sum(EJ) + sum(Ex),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy = ',EJ(1) + Ex(1),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy = ',EJ(2),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy = ',EJ(3) + Ex(2),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Hartree         energy = ',sum(EJ(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Hartree         energy = ',sum(EJ),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Hartree      aa energy = ',EJ(1),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Hartree      ab energy = ',EJ(2),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Hartree      bb energy = ',EJ(3),' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy = ',sum(Ex(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy = ',sum(Ex),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Exchange     a  energy = ',Ex(1),' au'
   write(*,'(A40,1X,F16.10,A3)') ' Exchange     b  energy = ',Ex(2),' au'
   write(*,'(A60)')              '---------------------------------------------'
diff --git a/src/LR/phULR_A.f90 b/src/LR/phULR_A.f90
index e4ce011..19fed4e 100644
--- a/src/LR/phULR_A.f90
+++ b/src/LR/phULR_A.f90
@@ -39,9 +39,9 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
   delta_dRPA = 0d0
   if(dRPA) delta_dRPA = 1d0
 
-!-----------------------------------------------
-! Build A matrix for spin-conserving transitions
-!-----------------------------------------------
+!----------------------------------------------
+! Build A matrix for spin-conserved transitions
+!----------------------------------------------
 
   if(ispin == 1) then 
 
@@ -59,10 +59,10 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
             Aph(ia,jb) = (eHF(a,1) - eHF(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
                        + lambda*ERI_aaaa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,b,j,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
     ! aabb block
 
@@ -77,10 +77,10 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
  
             Aph(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j) 
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
     ! bbaa block
 
@@ -95,10 +95,10 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
  
             Aph(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
     ! bbbb block
 
@@ -114,16 +114,16 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
             Aph(nSa+ia,nSa+jb) = (eHF(a,2) - eHF(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
                                + lambda*ERI_bbbb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,b,j,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
   end if
 
-!-----------------------------------------------
+!-----------------------------------------
 ! Build A matrix for spin-flip transitions
-!-----------------------------------------------
+!-----------------------------------------
 
   if(ispin == 2) then 
 
@@ -141,10 +141,10 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
             jb = jb + 1
             Aph(ia,jb) = (eHF(a,2) - eHF(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
                        - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,b,j,a)
-          end  do
-        end  do
-      end  do
-    end  do
+           end do
+         end do
+       end do
+     end do
 
     ! baba block
 
@@ -160,12 +160,11 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
             Aph(nSa+ia,nSa+jb) = (eHF(a,1) - eHF(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
                                - (1d0 - delta_dRPA)*lambda*ERI_aabb(b,i,a,j)
 
-          end  do
-        end  do
-      end  do
-    end  do
+           end do
+         end do
+       end do
+     end do
 
   end if
 
-
 end subroutine 
diff --git a/src/LR/phULR_B.f90 b/src/LR/phULR_B.f90
index a7bc6a9..a9a231e 100644
--- a/src/LR/phULR_B.f90
+++ b/src/LR/phULR_B.f90
@@ -57,10 +57,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
  
             Bph(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+           end do
+         end do
+       end do
+     end do
 
     ! aabb block
 
@@ -75,10 +75,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
  
             Bph(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b) 
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
     ! bbaa block
 
@@ -93,10 +93,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
  
             Bph(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
     ! bbbb block
 
@@ -111,10 +111,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
  
             Bph(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+         end do
+       end do
+     end do
+   end do
 
   end if
 
@@ -139,10 +139,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
  
             Bph(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,j,b,a)
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
 
     ! baab block
 
@@ -157,10 +157,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
  
             Bph(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(j,i,a,b)
 
-          end  do
-        end  do
-      end  do
-    end  do
+          end do
+        end do
+      end do
+    end do
   
   end if
 
diff --git a/src/LR/print_excitation_energies.f90 b/src/LR/print_excitation_energies.f90
index 2b43b21..4c67600 100644
--- a/src/LR/print_excitation_energies.f90
+++ b/src/LR/print_excitation_energies.f90
@@ -1,4 +1,4 @@
-subroutine print_excitation_energies(method,ispin,nS,Om)
+subroutine print_excitation_energies(method,manifold,nS,Om)
 
 ! Print excitation energies for a given spin manifold
 
@@ -8,36 +8,27 @@ subroutine print_excitation_energies(method,ispin,nS,Om)
 ! Input variables
 
   character(len=*),intent(in)        :: method
-  integer,intent(in)                 :: ispin
+  character(len=*),intent(in)        :: manifold
   integer,intent(in)                 :: nS
   double precision,intent(in)        :: Om(nS)
 
 ! Local variables
 
-  character(len=20)                  :: spin_manifold
   integer,parameter                  :: maxS = 20
-  integer                            :: ia
-
-  if(ispin == 1) spin_manifold = 'singlet'
-  if(ispin == 2) spin_manifold = 'triplet'
-  if(ispin == 3) spin_manifold = 'alpha-beta'
-  if(ispin == 4) spin_manifold = 'alpha-alpha' 
-  if(ispin == 5) spin_manifold = 'spin-conserved'
-  if(ispin == 6) spin_manifold = 'spin-flip'
-  if(ispin == 7) spin_manifold = 'beta-beta'
+  integer                            :: m
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------'
-  write(*,'(1X,A15,A15,A15,A9)') trim(method),' calculation: ',trim(spin_manifold),' manifold'
+  write(*,'(1X,A15,A15,A15,A9)') trim(method),' calculation: ',trim(manifold),' manifold'
   write(*,*)'-------------------------------------------------------------'
   write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
             '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
   write(*,*)'-------------------------------------------------------------'
 
-  do ia=1,min(maxS,nS)
+  do m=1,min(maxS,nS)
     write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & 
-      '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
-  enddo
+      '|',m,'|',Om(m),'|',Om(m)*HaToeV,'|'
+  end do
 
   write(*,*)'-------------------------------------------------------------'
   write(*,*)
diff --git a/src/RPA/crGRPA.f90 b/src/RPA/crGRPA.f90
index 3314b29..2f0a481 100644
--- a/src/RPA/crGRPA.f90
+++ b/src/RPA/crGRPA.f90
@@ -65,7 +65,7 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
   call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  call print_excitation_energies('crRPA@GHF',ispin,nS,Om)
+  call print_excitation_energies('crRPA@GHF','spinorbital',nS,Om)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   write(*,*)
diff --git a/src/RPA/crRRPA.f90 b/src/RPA/crRRPA.f90
index 31bbf53..2d7887e 100644
--- a/src/RPA/crRRPA.f90
+++ b/src/RPA/crRRPA.f90
@@ -73,7 +73,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
     if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
     call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    call print_excitation_energies('crRPA@RHF',ispin,nS,Om)
+    call print_excitation_energies('crRPA@RHF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   endif
@@ -88,7 +88,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
     if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
     call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    call print_excitation_energies('crRPA@HF',ispin,nS,Om)
+    call print_excitation_energies('crRPA@HF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   endif
diff --git a/src/RPA/phGRPA.f90 b/src/RPA/phGRPA.f90
index 7aae50f..21d6716 100644
--- a/src/RPA/phGRPA.f90
+++ b/src/RPA/phGRPA.f90
@@ -65,7 +65,7 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  call print_excitation_energies('phRPA@GHF',ispin,nS,Om)
+  call print_excitation_energies('phRPA@GHF','spinorbital',nS,Om)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   write(*,*)
diff --git a/src/RPA/phGRPAx.f90 b/src/RPA/phGRPAx.f90
index 9f7f5d7..1137322 100644
--- a/src/RPA/phGRPAx.f90
+++ b/src/RPA/phGRPAx.f90
@@ -65,7 +65,7 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  call print_excitation_energies('phRPAx@GHF',ispin,nS,Om)
+  call print_excitation_energies('phRPAx@GHF','spinorbital',nS,Om)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   write(*,*)
diff --git a/src/RPA/phRRPA.f90 b/src/RPA/phRRPA.f90
index b69e7fb..403e4f3 100644
--- a/src/RPA/phRRPA.f90
+++ b/src/RPA/phRRPA.f90
@@ -73,7 +73,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
     if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
     call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    call print_excitation_energies('phRPA@RHF',ispin,nS,Om)
+    call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   endif
@@ -88,7 +88,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
     if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
     call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    call print_excitation_energies('phRPA@RHF',ispin,nS,Om)
+    call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   endif
diff --git a/src/RPA/phRRPAx.f90 b/src/RPA/phRRPAx.f90
index 48415be..5021fd1 100644
--- a/src/RPA/phRRPAx.f90
+++ b/src/RPA/phRRPAx.f90
@@ -74,7 +74,7 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
     if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
     call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    call print_excitation_energies('phRPAx@HF',ispin,nS,Om)
+    call print_excitation_energies('phRPAx@HF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   endif
@@ -89,7 +89,7 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
     if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
   
     call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    call print_excitation_energies('phRPAx@HF',ispin,nS,Om)
+    call print_excitation_energies('phRPAx@HF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   endif
diff --git a/src/RPA/phURPA.f90 b/src/RPA/phURPA.f90
index d4af923..d5624d0 100644
--- a/src/RPA/phURPA.f90
+++ b/src/RPA/phURPA.f90
@@ -50,7 +50,7 @@ subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nB
   double precision              :: EcRPA(nspin)
 
 ! Hello world
-
+ 
   write(*,*)
   write(*,*)'***********************************'
   write(*,*)'* Unrestricted ph-RPA Calculation *'
@@ -87,7 +87,7 @@ subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nB
     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)
-    call print_excitation_energies('phRPA@UHF',5,nSt,Om)
+    call print_excitation_energies('phRPA@UHF','spin-flip',nSt,Om)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
@@ -112,7 +112,7 @@ subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nB
     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)
-    call print_excitation_energies('phRPA@UHF',6,nSt,Om)
+    call print_excitation_energies('phRPA@UHF','spin-conserved',nSt,Om)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
diff --git a/src/RPA/phURPAx.f90 b/src/RPA/phURPAx.f90
index 3b5017f..8dffc71 100644
--- a/src/RPA/phURPAx.f90
+++ b/src/RPA/phURPAx.f90
@@ -87,7 +87,7 @@ subroutine phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,n
     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)
-    call print_excitation_energies('phRPA@UHF',5,nSt,Om)
+    call print_excitation_energies('phRPA@UHF','spin-conserved',nSt,Om)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
@@ -112,7 +112,7 @@ subroutine phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,n
     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)
-    call print_excitation_energies('phRPA@UHF',6,nSt,Om)
+    call print_excitation_energies('phRPA@UHF','spin-flip',nSt,Om)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
diff --git a/src/RPA/ppGRPA.f90 b/src/RPA/ppGRPA.f90
index 132a10b..6a9c62a 100644
--- a/src/RPA/ppGRPA.f90
+++ b/src/RPA/ppGRPA.f90
@@ -66,8 +66,8 @@ subroutine ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,ENuc,EGHF,ERI,dipole_int,eHF)
 
 !   call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-  call print_excitation_energies('ppRPA@GHF (N+2)',ispin,nVV,Om1)
-  call print_excitation_energies('ppRPA@GHF (N-2)',ispin,nOO,Om2)
+  call print_excitation_energies('ppRPA@GHF','2p (spinorbital)',nVV,Om1)
+  call print_excitation_energies('ppRPA@GHF','2h (spinorbital)',nOO,Om2)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
diff --git a/src/RPA/ppRRPA.f90 b/src/RPA/ppRRPA.f90
index db5ee38..fc5dd77 100644
--- a/src/RPA/ppRRPA.f90
+++ b/src/RPA/ppRRPA.f90
@@ -78,8 +78,8 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E
 
 !   call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-    call print_excitation_energies('ppRPA@HF (N+2)',ispin,nVV,Om1)
-    call print_excitation_energies('ppRPA@HF (N-2)',ispin,nOO,Om2)
+    call print_excitation_energies('ppRPA@RHF','2p (singlet)',nVV,Om1)
+    call print_excitation_energies('ppRPA@RHF','2h (singlet)',nOO,Om2)
 
     deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp)
 
@@ -110,8 +110,8 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E
 
 !   call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-    call print_excitation_energies('ppRPA@HF (N+2)',ispin,nVV,Om1)
-    call print_excitation_energies('ppRPA@HF (N-2)',ispin,nOO,Om2)
+    call print_excitation_energies('ppRPA@RHF','2p (triplet)',nVV,Om1)
+    call print_excitation_energies('ppRPA@RHF','2h (triplet)',nOO,Om2)
 
     deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp)
 
diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90
index e652bf4..0c639c4 100644
--- a/src/RPA/ppURPA.f90
+++ b/src/RPA/ppURPA.f90
@@ -74,8 +74,8 @@ subroutine ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,E
              ERI_aabb,ERI_bbbb,Om1sc,X1sc,Y1sc, &
              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)
+  call print_excitation_energies('ppRPA@UHF','2p (alpha-beta)',nP_sc,Om1sc)
+  call print_excitation_energies('ppRPA@UHF','2h (alpha-beta)',nH_sc,Om2sc)
 
 !alpha-alpha block
 
@@ -100,8 +100,8 @@ subroutine ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,E
              ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf, &
              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)
+  call print_excitation_energies('ppRPA@UHF','2h (alpha-alpha)',nP_sf,Om1sf)
+  call print_excitation_energies('ppRPA@UHF','2p (alpha-alpha)',nH_sf,Om2sf)
 
   deallocate(Om1sf,X1sf,Y1sf,Om2sf,X2sf,Y2sf)
 
@@ -120,8 +120,8 @@ subroutine ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,E
              ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf,&
              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)
+  call print_excitation_energies('ppRPA@UHF','2p (beta-beta)',nP_sf,Om1sf)
+  call print_excitation_energies('ppRPA@UHF','2h (beta-beta)',nH_sf,Om2sf)
 
   EcRPA(2) = 3d0*EcRPA(2)