From 4693d18fc03042870347feb7047338e99a84efcd Mon Sep 17 00:00:00 2001
From: "Dr. Mathieu Taillefumier" <mathieu.taillefumier@free.fr>
Date: Thu, 24 Nov 2022 09:59:12 +0100
Subject: [PATCH] Add mixing factor 2 in Hubbard U one-electron contribution

---
 src/dft/energy.cpp                       | 9 ++++-----
 src/hubbard/hubbard_potential_energy.cpp | 7 +++++++
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/dft/energy.cpp b/src/dft/energy.cpp
index dde7271a0..2156abe6f 100644
--- a/src/dft/energy.cpp
+++ b/src/dft/energy.cpp
@@ -165,7 +165,7 @@ 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;
@@ -173,8 +173,7 @@ total_energy(Simulation_context const& ctx, K_point_set const& kset, Density con
     }
 
     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;
 }
diff --git a/src/hubbard/hubbard_potential_energy.cpp b/src/hubbard/hubbard_potential_energy.cpp
index 61045968e..0d01b5c5f 100644
--- a/src/hubbard/hubbard_potential_energy.cpp
+++ b/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;