Add mixing factor 2 in Hubbard U one-electron contribution

src/dft/energy.cpp

@@ -165,16 +165,15 @@ total_energy(Simulation_context const& ctx, K_point_set const& kset, Density con
 
         case electronic_structure_method_t::pseudopotential: {
             tot_en = (kset.valence_eval_sum() - energy_vxc(density, potential) - energy_bxc(density, potential) -
-                      potential.PAW_one_elec_energy(density)) -
+                      potential.PAW_one_elec_energy(density) - one_electron_energy_hubbard(density, potential)) -
                      0.5 * energy_vha(potential) + energy_exc(density, potential) + potential.PAW_total_energy() +
                      ewald_energy + kset.entropy_sum();
             break;
         }
     }
 
     if (ctx.hubbard_correction()) {
-        tot_en += ::sirius::energy(density.occupation_matrix());
-        tot_en -= ::sirius::one_electron_energy_hubbard(density, potential);
+        tot_en += ::sirius::hubbard_energy(density);
     }
 
     return tot_en;
@@ -194,7 +193,7 @@ double
 one_electron_energy(Density const& density, Potential const& potential)
 {
     return energy_vha(potential) + energy_vxc(density, potential) + energy_bxc(density, potential) +
-           potential.PAW_one_elec_energy(density) + one_electron_energy_hubbard(density, potential);
+      potential.PAW_one_elec_energy(density) + one_electron_energy_hubbard(density, potential);
 }
 
 double
@@ -213,7 +212,7 @@ energy_potential(Density const& density, Potential const& potential)
     double e =
         energy_veff(density, potential) + energy_bxc(density, potential) + potential.PAW_one_elec_energy(density);
     if (potential.ctx().hubbard_correction()) {
-        e += ::sirius::energy(density.occupation_matrix());
+        e += ::sirius::hubbard_energy(density);
     }
     return e;
 }

src/hubbard/hubbard_potential_energy.cpp

@@ -188,6 +188,9 @@ calculate_energy_collinear_nonlocal(Simulation_context const& ctx__, const int i
         }
     }
 
+    if (ctx__.num_spins() == 1) {
+        hubbard_energy *= 2.0;
+    }
     return -0.5 * hubbard_energy;
 }
 
@@ -631,6 +634,10 @@ one_electron_energy_hubbard(Hubbard_matrix const& om__, Hubbard_matrix const& pm
               }
             }
         }
+
+        if (ctx.num_spins() == 1) {
+          tmp *= 2.0;
+        }
         return std::real(tmp);
     }
     return 0.0;