Merge branch 'develop' into develop

Author: ESROAMER

source/source_esolver/esolver_ks.cpp

@@ -13,7 +13,7 @@
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_io/module_output/output_log.h" // use write_head
 #include "source_estate/elecstate_print.h" // print_etot
-#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07
+#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 2025-11-07
 #include "source_hamilt/module_xc/general_exx_info.h" // for init_general_exx_info
 
 namespace ModuleESolver

source/source_esolver/esolver_ks.h

@@ -8,7 +8,7 @@
 #include "source_hamilt/hamilt.h"                      // use Hamiltonian
 #include "source_hamilt/hamilt_base.h"                 // use Hamiltonian base class
 #include "source_hamilt/module_xc/general_exx_info.h"  // ESolver owns General_Exx_Info value
-#include "source_lcao/module_dftu/dftu.h"              // mohan add 20251107
+#include "source_lcao/module_dftu/dftu_lcao.h"              // mohan add 20251107
 #include "source_pw/module_pwdft/vnl_pw.h"
 
 namespace ModuleESolver

source/source_esolver/esolver_ks_lcao.cpp

@@ -508,6 +508,14 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     // 1) calculate the local occupation number matrix and energy correction in DFT+U
     finish_dftu_lcao<TK>(iter, conv_esolver, this->inp_->dft_plus_u, this->inp_->out_chg[0], &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, PARAM.globalv.global_out_dir, this->inp_->nspin, PARAM.globalv.npol);
 
+    // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState.
+    // Covers both dft_plus_u==1 (new method, energy accumulated by DFTU::contributeHR
+    // via cal_v_of_u) and dft_plus_u==2 (old method, energy from cal_energy_correction).
+    if (this->inp_->dft_plus_u)
+    {
+        this->pelec->set_dftu_energy(this->dftu.get_energy());
+    }
+
     // 2) for deepks, calculate delta_e, output labels during electronic steps
     this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, *this->inp_);
 

source/source_esolver/esolver_ks_pw.cpp

@@ -210,6 +210,12 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
                           ucell,
                           this->p_chgmix,
                           this->kv.isk.data());
+
+    // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState
+    if (this->inp_->dft_plus_u)
+    {
+        this->pelec->set_dftu_energy(this->dftu.get_energy());
+    }
 }
 
 // Temporary, it should be replaced by hsolver later.

source/source_esolver/lcao_others.cpp

@@ -2,7 +2,7 @@
 #include "source_cell/cal_ux.h"
 #include "source_estate/module_charge/symm_rho.h"
 #include "source_lcao/hamilt_lcao.h"
-#include "source_lcao/module_dftu/dftu.h"
+#include "source_lcao/module_dftu/dftu_lcao.h"
 #include "source_hamilt/module_gint/gint.h"
 #include "source_base/formatter.h"
 #include "source_base/timer.h"

source/source_estate/elecstate.h

@@ -131,6 +131,7 @@ class ElecState
     }
 
     double get_dftu_energy();
+    void set_dftu_energy(double e) { dftu_energy_ = e; }
     double get_local_pp_energy();
 
     fenergy f_en; ///< energies contribute to the total free energy
@@ -148,6 +149,9 @@ class ElecState
   public:
 
     bool skip_weights = false;
+
+  private:
+    double dftu_energy_ = 0.0; ///< DFT+U energy, set by ESolver via set_dftu_energy()
 };
 
 /**

source/source_estate/estate_e_terms.cpp

@@ -3,7 +3,6 @@
 #include "source_estate/module_pot/h_hartree_pw.h"
 #include "source_estate/module_pot/efield.h"
 #include "source_estate/module_pot/gatefield.h"
-#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06
 
 namespace elecstate
 {
@@ -35,7 +34,7 @@ double ElecState::get_solvent_model_Acav()
 
 double ElecState::get_dftu_energy()
 {
-    return Plus_U::get_energy();
+    return dftu_energy_;
 }
 
 double ElecState::get_local_pp_energy()

source/source_estate/module_charge/chgmixing.cpp

@@ -1,6 +1,6 @@
 #include "source_estate/module_charge/chgmixing.h"
 #include "source_estate/update_pot.h"
-#include "source_lcao/module_dftu/dftu.h"
+#include "source_lcao/module_dftu/dftu_lcao.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
 
 void module_charge::chgmixing_ks(const int iter, // scf iteration number
@@ -119,7 +119,7 @@ void module_charge::chgmixing_ks(const int iter, // scf iteration number
 
 void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
         Charge_Mixing* p_chgmix, // charge mixing class
-        Plus_U &dftu, // mohan add 2025-11-06
+        Plus_U_Base& dftu,
 		const Input_para& inp) // input parameters
 {
     ModuleBase::TITLE("module_charge", "chgmixing_ks_pw");
@@ -167,7 +167,7 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
 					std::cout << " U-Ramping! Current U = ";
 					for (int i = 0; i < dftu.get_num_u_types(); i++)
 					{
-						std::cout << dftu.get_hubbard_u(i) * ModuleBase::Ry_to_eV << " ";
+						std::cout << dftu.get_u_current(i) * ModuleBase::Ry_to_eV << " ";
 					}
 					std::cout << " eV " << std::endl;
 				}
@@ -180,7 +180,7 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
 
 void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number
         Charge_Mixing* p_chgmix, // charge mixing class
-        Plus_U &dftu, // mohan add 2025-11-06
+        Plus_U_Base& dftu,
         const int nnr, // dimension of density matrix
 		const Input_para& inp) // input parameters
 {
@@ -202,7 +202,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number
             std::cout << " U-Ramping! Current U = ";
             for (int i = 0; i < dftu.get_num_u_types(); i++)
             {
-                std::cout << dftu.get_hubbard_u(i) * ModuleBase::Ry_to_eV << " ";
+                std::cout << dftu.get_u_current(i) * ModuleBase::Ry_to_eV << " ";
             }
             std::cout << " eV " << std::endl;
         }
@@ -221,7 +221,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number
                 std::cout << " U-Ramping! Current U = ";
                 for (int i = 0; i < dftu.get_num_u_types(); i++)
                 {
-                    std::cout << dftu.get_hubbard_u(i) * ModuleBase::Ry_to_eV << " ";
+                    std::cout << dftu.get_u_current(i) * ModuleBase::Ry_to_eV << " ";
                 }
                 std::cout << " eV " << std::endl;
             }

source/source_estate/module_charge/chgmixing.h

@@ -6,7 +6,7 @@
 #include "source_estate/module_charge/charge_mixing.h" // use p_chgmix
 #include "source_io/module_parameter/input_parameter.h" // use Input_para
 #include "source_cell/unitcell.h"
-#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06
+#include "source_pw/module_pwdft/dftu_base.h"
 
 namespace module_charge
 {
@@ -28,12 +28,12 @@ void chgmixing_ks(const int iter, // scf iteration number
 
 void chgmixing_ks_pw(const int iter,
         Charge_Mixing* p_chgmix,
-        Plus_U &dftu, // mohan add 2025-11-06
+        Plus_U_Base& dftu,
 		const Input_para& inp); // input parameters
 
 void chgmixing_ks_lcao(const int iter, // scf iteration number
         Charge_Mixing* p_chgmix, // charge mixing class
-        Plus_U &dftu, // mohan add 2025-11-06
+        Plus_U_Base& dftu,
         const int nnr, // dimension of density matrix
 		const Input_para& inp); // input parameters
 

source/source_hsolver/test/diago_mock.h

@@ -577,7 +577,7 @@ template<> hamilt::HamiltPW<double>::HamiltPW(
 		ModulePW::PW_Basis_K* wfc_basis,
 		K_Vectors* pkv,
 		pseudopot_cell_vnl* ppcell,
-		Plus_U* p_dftu, // mohan add 20251108
+		Plus_U_Base* p_dftu,
 		const UnitCell* ucell,
 		const General_Exx_Info* exx_info)
 {
@@ -599,7 +599,7 @@ template<> hamilt::HamiltPW<std::complex<double>>::HamiltPW(
 		ModulePW::PW_Basis_K* wfc_basis,
 		K_Vectors* pkv,
 		pseudopot_cell_vnl* ppcell,
-		Plus_U* p_dftu, // mohan add 20251108
+		Plus_U_Base* p_dftu,
 		const UnitCell* ucell,
 		const General_Exx_Info* exx_info)
 {
@@ -621,7 +621,7 @@ template<> hamilt::HamiltPW<std::complex<float>>::HamiltPW(
 		ModulePW::PW_Basis_K* wfc_basis,
 		K_Vectors* pkv,
 		pseudopot_cell_vnl* ppcell,
-		Plus_U* p_dftu, // mohan add 20251108
+		Plus_U_Base* p_dftu,
 		const UnitCell* ucell,
 		const General_Exx_Info* exx_info)
 {