Refactor DFT+U: rename and de-static Plus_U members, step 3 (deepmodeling#7859)

* Refactor DFT+U: rename and de-static Plus_U members

- Rename U -> u_current, U0 -> u_target, omc -> occ_mat_ctrl, Yukawa -> use_yukawa
- Rename accessors: get_hubbard_u -> get_u_current, get_hubbard_u0 -> get_u_target
- Remove static qualifier from all 9 Plus_U_Base members (u_current, u_target,
  occ_mat_ctrl, use_yukawa, energy_u, orbital_corr, uramping, mixing_dftu, nspin)
- Convert static accessors (get_energy, set_energy, get_u_current, etc.) and
  is_mixing_enabled/enable_mixing to non-static member functions
- Extend Plus_U_Base::init_base() and Plus_U::init() signatures to accept
  hubbard_u, uramping, occ_mat_ctrl, mixing_dftu; move initialization logic
  out of input_conv.cpp into init_base() so members are set after the Plus_U
  instance is constructed (no longer relies on static state)
- Remove redundant Plus_U assignments from input_conv.cpp; drop its include
  of dftu.h (dependencies made explicit at call sites)
- Introduce ElecState::set_dftu_energy() / dftu_energy_ member so
  get_dftu_energy() no longer depends on Plus_U::get_energy() static state;
  ESolver_KS_PW and ESolver_KS_LAO push the energy after DFT+U update
- Replace PARAM.inp.* with inp.* in setup_pot.cpp (inp is already a parameter)
- Rename Input_para::omc -> occ_mat_ctrl (INPUT keyword 'omc' preserved for
  backward compatibility)

Verification: pending (compilation and tests to be run by user)

* fix bug

* fix bugs

* add an example to test the old DFT+U method with only scf energy, not force and stress

* update CASES_CPU.txt

* delete useless members

* Move LCAO-only DFT+U members from Plus_U_Base to Plus_U

Seven members (yukawa_lambda, npol, nlocal, gamma_only_local, ks_solver, cal_force, cal_stress) were only read by the LCAO subclass (Plus_U) and its component files under source_lcao/module_dftu/. They are now declared in Plus_U and assigned in Plus_U::init, so Plus_U_Base no longer carries LCAO-specific mutable workflow state.

Accordingly shorten Plus_U_Base::init_base signature by six parameters (yukawa_lambda, nlocal, gamma_only_local, ks_solver, cal_force, cal_stress); update the two call sites (Plus_U::init in dftu.cpp, PW path in setup_pot.cpp). Plus_U::init keeps its original signature so external callers (lcao_set.cpp) are unaffected.

* Inline DFTU template header into dftu_lcao.h

dftu.hpp only provided the empty primary template for hamilt::DFTU. Move it into dftu_lcao.h (the sole includer), guarded by __DFTUTEMPLATE, matching the convention used by sibling LCAO operators (nonlocal.h, ekinetic.h, meta_lcao.h, veff_lcao.h, etc.) that inline the primary template instead of using a separate .hpp file. This also aligns with AGENTS.md rule 4 (avoid .hpp implementation headers).

* Rename DFT+U LCAO files to reflect roles

Three LCAO DFT+U files had confusing names that hid which was the main implementation, the new operator, and the legacy operator. Rename for clarity:

  - dftu.h/.cpp -> dftu_lcao.h/.cpp   : Plus_U class (LCAO DFT+U physics, subclass of Plus_U_Base in dftu_base.h). The name now mirrors dftu_base.h.

  - dftu_lcao.h/.cpp -> dftu_lcao_op.h/.cpp : hamilt::DFTU operator (new method, dft_plus_u==1, onsite_radius projections). Operators now carry the _op suffix.

  - op_dftu_lcao.h/.cpp -> dftu_lcao_op_legacy.h/.cpp : hamilt::OperatorDFTU operator (old method, dft_plus_u==2, NAO first zeta projections, test-only). Renamed to mark it as legacy and group it with the new operator.

Class names are unchanged in this step (Plus_U, hamilt::DFTU, hamilt::OperatorDFTU); only filenames and #include paths change. Updated all ~40 includers across source_lcao, source_pw, source_esolver, source_estate, source_io, plus the module and test CMakeLists (the dftu_lcao_test links the new-operator cpp ../dftu_lcao_op.cpp, not the Plus_U cpp).

* refactor: eliminate reverse dependency of dftu_lcao.h from PW/base layers

Replace Plus_U (LCAO-derived class) with Plus_U_Base (PW-layer base class) in
all non-LCAO interfaces. All DFT+U calls in PW/estate/test layers only use
base-class members, so the derived-type reference was unnecessary.

Changes:
- Remove #include "source_lcao/module_dftu/dftu_lcao.h" from PW/base headers
- Use #include "source_pw/module_pwdft/dftu_base.h" instead
- Replace Plus_U&/Plus_U* with Plus_U_Base&/Plus_U_Base* in:
  * source_pw/module_pwdft: force_pw, hamilt_pw, onsite_proj, op_pw_proj,
    stress_func, stress_onsite, stress_pw (6 headers + 7 impls)
  * source_estate/module_charge: chgmixing (PW + LCAO mixing functions)
  * source_pw/module_stodft/test: test_hamilt_sto, test_sto_tool
  * source_hsolver/test: diago_mock.h (3 explicit specializations)

Plus_U (derived) instances held by callers (e.g. ESolver_KS::dftu) continue
to work via implicit pointer/reference upcast to Plus_U_Base. No behavioral
change.

---------

Co-authored-by: abacus_fixer <mohanchen@pku.eud.cn>

Author: mohanchen

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)
 {