diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90
index f27e907..56279ff 100644
--- a/src/eDFT/print_unrestricted_individual_energy.f90
+++ b/src/eDFT/print_unrestricted_individual_energy.f90
@@ -1,4 +1,5 @@
-subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD)
+subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, & 
+                                                Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD)
 
 ! Print individual energies for eDFT calculation
 
@@ -7,26 +8,29 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux
 
 ! Input variables
 
-  integer,intent(in)                 :: nEns
-  double precision,intent(in)        :: ENuc
-  double precision,intent(in)        :: Ew
-  double precision,intent(in)        :: ET(nspin,nEns)
-  double precision,intent(in)        :: EV(nspin,nEns)
-  double precision,intent(in)        :: EH(nsp,nEns)
-  double precision,intent(in)        :: Ex(nspin,nEns)
-  double precision,intent(in)        :: Ec(nsp,nEns)
-  double precision,intent(in)        :: Eaux(nspin,nEns)
-  double precision,intent(in)        :: ExDD(nspin,nEns)
-  double precision,intent(in)        :: EcDD(nsp,nEns)
-  double precision,intent(in)        :: E(nEns)
+  integer,intent(in)            :: nEns
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: Ew
+  double precision,intent(in)   :: ET(nspin,nEns)
+  double precision,intent(in)   :: EV(nspin,nEns)
+  double precision,intent(in)   :: EH(nsp,nEns)
+  double precision,intent(in)   :: Ex(nspin,nEns)
+  double precision,intent(in)   :: Ec(nsp,nEns)
+  double precision,intent(in)   :: Eaux(nspin,nEns)
+  double precision              :: LZH(nsp)
+  double precision              :: LZx(nspin)
+  double precision              :: LZc(nsp)
+  double precision,intent(in)   :: ExDD(nspin,nEns)
+  double precision,intent(in)   :: EcDD(nsp,nEns)
+  double precision,intent(in)   :: E(nEns)
 
-  double precision,intent(in)        :: OmH(nEns)
-  double precision,intent(in)        :: Omx(nEns)
-  double precision,intent(in)        :: Omc(nEns)
-  double precision,intent(in)        :: Omaux(nEns)
-  double precision,intent(in)        :: OmxDD(nEns)
-  double precision,intent(in)        :: OmcDD(nEns)
-  double precision,intent(in)        :: Om(nEns)
+  double precision,intent(in)   :: OmH(nEns)
+  double precision,intent(in)   :: Omx(nEns)
+  double precision,intent(in)   :: Omc(nEns)
+  double precision,intent(in)   :: Omaux(nEns)
+  double precision,intent(in)   :: OmxDD(nEns)
+  double precision,intent(in)   :: OmcDD(nEns)
+  double precision,intent(in)   :: Om(nEns)
 
 ! Local variables
 
@@ -134,6 +138,26 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux
   write(*,'(A60)')           '-------------------------------------------------'
   write(*,*)
 
+!------------------------------------------------------------------------
+! Print Levy-Zahariev shifts
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A60)')           ' LEVY-ZAHARIEV SHIFTS CONTRIBUTIONS'
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+  write(*,'(A40,F16.10,A3)') '   H Levy-Zahariev shift: ',sum(LZH(:)),' au'
+  write(*,'(A40,F16.10,A3)') '   x Levy-Zahariev shift: ',sum(LZx(:)),' au'
+  write(*,'(A40,F16.10,A3)') '   c Levy-Zahariev shift: ',sum(LZc(:)),' au'
+  write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZx(:)),' au'
+  write(*,*)
+  write(*,'(A40,F16.10,A3)') '   H Levy-Zahariev shift: ',sum(LZH(:))*HaToeV,' eV'
+  write(*,'(A40,F16.10,A3)') '   x Levy-Zahariev shift: ',sum(LZx(:))*HaToeV,' eV'
+  write(*,'(A40,F16.10,A3)') '   c Levy-Zahariev shift: ',sum(LZc(:))*HaToeV,' eV'
+  write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZx(:)))*HaToeV,' eV'
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+
 !------------------------------------------------------------------------
 ! Compute derivative discontinuities
 !------------------------------------------------------------------------
diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90
index 0227c75..7846bff 100644
--- a/src/eDFT/unrestricted_individual_energy.f90
+++ b/src/eDFT/unrestricted_individual_energy.f90
@@ -232,6 +232,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 ! Dump results
 !------------------------------------------------------------------------
 
-  call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD)
+  call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, & 
+                                            Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD)
 
 end subroutine unrestricted_individual_energy