diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90
index 7eb6221..c91d9fc 100644
--- a/src/eDFT/US51_lda_exchange_individual_energy.f90
+++ b/src/eDFT/US51_lda_exchange_individual_energy.f90
@@ -48,7 +48,7 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex
 
             rI = max(0d0,rho(iG,ispin,iEns))
 
-            if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(e+dedr*r)*rI
+            if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(e + dedr*r)*rI
 
           end do
 
diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90
index 05447ea..99b3159 100644
--- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90
+++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90
@@ -144,7 +144,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZ
           if(raI > threshold) then 
 
             Ec(1,iEns) = Ec(1,iEns) + weight(iG)*(ec_z + decdra*ra)*raI
-            Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdra*rb*raI
+            if(rb > threshold) Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdra*rb*raI
 
           end if
 
@@ -152,6 +152,10 @@ subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZ
 
       end if
 
+      ! up-down contribution
+
+      if(ra > threshold .and. rb > threshold) LZc(2) = LZc(2) -weight(iG)*(decdra + decdrb)*ra*rb
+
       ! spin-down contribution
 
       if(rb > threshold) then 
@@ -165,7 +169,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZ
           if(rbI > threshold) then 
 
             Ec(3,iEns) = Ec(3,iEns) + weight(iG)*(ec_z + decdrb*rb)*rbI
-            Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdrb*ra*rbI
+            if(ra > threshold) Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdrb*ra*rbI
 
           end if
 
diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90
index 71bf04e..283ea01 100644
--- a/src/eDFT/unrestricted_individual_energy.f90
+++ b/src/eDFT/unrestricted_individual_energy.f90
@@ -167,8 +167,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
               + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
     end do
 
-    print*,E
-
   else 
 
     do iEns=1,nEns
@@ -179,14 +177,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
   end if
   
-  do iEns=1,nEns
-    E(iEns)   = sum(ET(:,iEns)) + sum(EV(:,iEns))                   & 
-              + sum(EH(:,iEns)) + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) & 
-              + kappa(iEns)*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)))            &
-              + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
-  end do
+! do iEns=1,nEns
+!   E(iEns)   = sum(ET(:,iEns)) + sum(EV(:,iEns))                     & 
+!             + sum(EH(:,iEns)) + sum(Ex(:,iEns)) + sum(Ec(:,iEns))   & 
+!             + sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) &
+!             + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
+! end do
 
-  print*,E
+!   print*,E
 
 !------------------------------------------------------------------------
 ! Excitation energies
@@ -196,29 +194,29 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
     Om(iEns) = E(iEns) - E(1)
 
-    OmH(iEns)    = sum(EH(:,iEns)) - sum(EH(:,1)) 
-    Omx(iEns)    = sum(Ex(:,iEns)) - sum(Ex(:,1))
-    Omc(iEns)    = sum(Ec(:,iEns)) - sum(Ec(:,1))
+    OmH(iEns)   = sum(EH(:,iEns)) - sum(EH(:,1)) 
+    Omx(iEns)   = sum(Ex(:,iEns)) - sum(Ex(:,1))
+    Omc(iEns)   = sum(Ec(:,iEns)) - sum(Ec(:,1))
 
-    Omaux(iEns)  = sum(Eaux(:,iEns)) - sum(Eaux(:,1))
+    Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1))
 
-    OmxDD(iEns)  = sum(ExDD(:,iEns)) - sum(ExDD(:,1))
-    OmcDD(iEns)  = sum(EcDD(:,iEns)) - sum(EcDD(:,1))
+    OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1))
+    OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1))
 
   end do
 
   if(doNcentered) then
 
     do iEns=1,nEns
-      OmH(iEns)    = OmH(iEns) + (kappa(iEns) - kappa(1))*sum(LZH(:))
-      Omx(iEns)    = Omx(iEns) + (kappa(iEns) - kappa(1))*sum(LZx(:))
-      Omc(iEns)    = Omc(iEns) + (kappa(iEns) - kappa(1))*sum(LZc(:))
+      OmH(iEns)   = OmH(iEns) + (kappa(iEns) - kappa(1))*sum(LZH(:))
+      Omx(iEns)   = Omx(iEns) + (kappa(iEns) - kappa(1))*sum(LZx(:))
+      Omc(iEns)   = Omc(iEns) + (kappa(iEns) - kappa(1))*sum(LZc(:))
 
-      Omaux(iEns)  = Omaux(iEns) &
-                   + (kappa(iEns) - kappa(1))*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))) 
+      Omaux(iEns) = Omaux(iEns) &
+                  + (kappa(iEns) - kappa(1))*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))) 
 
-      OmxDD(iEns)  = kappa(iEns)*sum(ExDD(:,iEns)) - kappa(1)*sum(ExDD(:,1))
-      OmcDD(iEns)  = kappa(iEns)*sum(EcDD(:,iEns)) - kappa(1)*sum(EcDD(:,1))
+      OmxDD(iEns) = kappa(iEns)*sum(ExDD(:,iEns)) - kappa(1)*sum(ExDD(:,1))
+      OmcDD(iEns) = kappa(iEns)*sum(EcDD(:,iEns)) - kappa(1)*sum(EcDD(:,1))
   
     end do