Refactor DFT+U: rename and de-static Plus_U members, step 3

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

@@ -506,6 +506,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

@@ -209,6 +209,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)
 {

source/source_io/module_ctrl/ctrl_scf_lcao.h

@@ -9,7 +9,7 @@
 #include "source_estate/module_dm/density_matrix.h"    // mohan add 2025-11-04
 #include "source_hamilt/module_surchem/surchem.h"       // use surchem (for dH veff pots)
 #include "source_lcao/hamilt_lcao.h"                   // use hamilt::HamiltLCAO<TK, TR>
-#include "source_lcao/module_dftu/dftu.h"              // mohan add 20251107
+#include "source_lcao/module_dftu/dftu_lcao.h"              // mohan add 20251107
 #include "source_lcao/module_rdmft/rdmft.h"            // use RDMFT codes
 #include "source_lcao/setup_deepks.h"                  // for deepks, mohan add 20251008
 #include "source_lcao/setup_exx.h"                     // for exx, mohan add 20251008

source/source_io/module_hs/output_mat_sparse.h

@@ -6,7 +6,7 @@
 #include "source_cell/klist.h"
 #include "source_hamilt/hamilt.h"
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
-#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107
+#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107
 
 namespace ModuleIO
 {

source/source_io/module_hs/write_hs_r.h

@@ -6,7 +6,7 @@
 #include "source_cell/klist.h"
 #include "source_hamilt/hamilt.h"
 #include "source_lcao/lcao_hs_arrays.hpp"
-#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107
+#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107
 
 #ifdef __EXX
 #include "RI/global/Tensor.h" // for RI::Tensor

source/source_io/module_hs/write_vxc.hpp

@@ -4,7 +4,7 @@
 #include "source_base/parallel_reduce.h"
 #include "source_base/module_container/base/third_party/blas.h"
 #include "source_base/module_external/scalapack_connector.h"
-#include "source_lcao/module_dftu/op_dftu_lcao.h"
+#include "source_lcao/module_dftu/dftu_lcao_op_legacy.h"
 #include "source_lcao/module_operator_lcao/veff_lcao.h"
 #include "source_hamilt/module_xc/exx_info.h"
 #ifdef __EXX

source/source_io/module_hs/write_vxc_r.hpp

@@ -2,7 +2,7 @@
 #define __WRITE_VXC_R_H_
 #include "source_io/module_parameter/parameter.h"
 #include "source_io/module_hs/write_hs_sparse.h"
-#include "source_lcao/module_dftu/op_dftu_lcao.h"
+#include "source_lcao/module_dftu/dftu_lcao_op_legacy.h"
 #include "source_lcao/module_operator_lcao/veff_lcao.h"
 #include "source_lcao/spar_hsr.h"
 #ifdef __EXX

source/source_io/module_parameter/input_conv.cpp

@@ -20,7 +20,6 @@
 #include "source_lcao/module_ri/exx_abfs_jle.h"
 #endif
 
-#include "source_lcao/module_dftu/dftu.h"
 #ifdef __LCAO
 #include "source_basis/module_ao/orb_read.h"
 #include "source_lcao/force_stress_lcao.h"
@@ -93,19 +92,9 @@ void Input_Conv::Convert()
     // iteration (1/3)
     //----------------------------------------------------------
 
-    if (PARAM.inp.dft_plus_u)
-    {
-        Plus_U::omc = PARAM.inp.omc;
-        Plus_U::orbital_corr = PARAM.inp.orbital_corr;
-        Plus_U::uramping = PARAM.globalv.uramping;
-        Plus_U::mixing_dftu = PARAM.inp.mixing_dftu;
-        Plus_U::U = PARAM.globalv.hubbard_u;
-        Plus_U::U0 = PARAM.globalv.hubbard_u;
-        if (PARAM.globalv.uramping > 0.01)
-        {
-            ModuleBase::GlobalFunc::ZEROS(Plus_U::U.data(), PARAM.inp.ntype);
-        }
-    }
+    // Note: DFT+U static members (u_current, u_target, occ_mat_ctrl, etc.)
+    // are now initialized inside Plus_U_Base::init_base() which is called
+    // from setup_pot.cpp (PW) and lcao_set.cpp (LCAO). Mohan refactor 2025-11.
 
     //----------------------------------------------------------
     // Yu Liu add 2022-05-18

source/source_io/module_parameter/input_parameter.h

@@ -626,7 +626,7 @@ struct Input_para
     bool yukawa_potential = false;         ///< default: false
     double yukawa_lambda = -1.0;           ///< default: -1.0, which means we calculate lambda
     double uramping_eV = -1.0;             ///< U-Ramping method (eV)
-    int omc = 0;                           ///< the mode of occupation matrix control
+    int occ_mat_ctrl = 0;                  ///< the mode of occupation matrix control
     double onsite_radius = 0.0;            ///< radius of the sphere for onsite projection (Bohr)
     std::vector<double> hubbard_u_eV = {}; ///< Hubbard Coulomb interaction parameter U(ev)
     std::vector<int> orbital_corr = {};    ///< which correlated orbitals need corrected ; d:2 ,f:3, do not

source/source_io/module_parameter/read_inp_exx_dftu.cpp

@@ -818,7 +818,7 @@ void ReadInput::item_dftu()
 [NOTE] The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward.)";
         item.default_value = "0";
         item.unit = "";
-        read_sync_int(input.omc);
+        read_sync_int(input.occ_mat_ctrl);
         this->add_item(item);
     }
     {

source/source_io/test/for_testing_input_conv.h

@@ -13,7 +13,7 @@
 #include "source_io/module_restart/restart.h"
 #include "source_io/module_unk/berryphase.h"
 #include "source_lcao/force_stress_lcao.h"
-#include "source_lcao/module_dftu/dftu.h"
+#include "source_lcao/module_dftu/dftu_lcao.h"
 #include "source_md/md_func.h"
 #include "source_pw/module_pwdft/stru_fac.h"
 #include "source_pw/module_pwdft/vnl_pw.h"

source/source_io/test/read_input_ptest.cpp

@@ -358,7 +358,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_FALSE(param.inp.yukawa_potential);
     EXPECT_DOUBLE_EQ(param.inp.yukawa_lambda, -1.0);
     EXPECT_EQ(param.inp.onsite_radius, 0.0);
-    EXPECT_EQ(param.inp.omc, 0);
+    EXPECT_EQ(param.inp.occ_mat_ctrl, 0);
     EXPECT_FALSE(param.inp.dft_plus_dmft);
     EXPECT_FALSE(param.inp.rpa);
     EXPECT_FALSE(param.inp.rpa_out_vel);

source/source_lcao/force_stress_lcao.cpp

@@ -1,7 +1,7 @@
 #include "force_stress_lcao.h"
 
 #include "source_base/parallel_reduce.h"
-#include "source_lcao/module_dftu/dftu.h" //Quxin add for DFT+U on 20201029
+#include "source_lcao/module_dftu/dftu_lcao.h" //Quxin add for DFT+U on 20201029
 #include "source_io/module_output/output_log.h"
 #include "source_io/module_parameter/parameter.h"
 // new
@@ -20,7 +20,7 @@
 #include "source_lcao/module_deepks/lcao_deepks_io.h" // mohan add 2024-07-22
 #include "source_lcao/module_deepks/deepks_force.h"
 #endif
-#include "source_lcao/module_dftu/dftu_lcao.h"
+#include "source_lcao/module_dftu/dftu_lcao_op.h"
 #include "source_lcao/module_operator_lcao/dspin_lcao.h"
 #include "source_lcao/module_operator_lcao/nonlocal.h"
 #include "source_lcao/module_operator_lcao/ekinetic.h"
@@ -426,7 +426,32 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         }
         if (PARAM.inp.dft_plus_u == 2)
         {
-            // Old DFT+U implementation (dft_plus_u==2) still needs ForceStressArrays
+            // The legacy dft_plus_u==2 force/stress path is currently broken.
+            //
+            // Background: Plus_U::force_stress relies on ForceStressArrays
+            // members DSloc_x/y/z (gamma_only) or DSloc_Rx/Ry/Rz (multik)
+            // and DH_r being pre-allocated and filled with dS/dR data by the
+            // main force flow (formerly ForceLcaoGamma::ftable). The DFT+U
+            // step 2 refactor (commit 70c54c9d5a, 2026-01-23) removed the
+            // main-flow ForceStressArrays because the operator-based force
+            // calculation no longer needs it, but the legacy dft_plus_u==2
+            // path still depends on it. The local fsr_dftu below is declared
+            // without allocating those arrays, so any call into
+            // cal_force_gamma / cal_stress_gamma / folding_matrix_k would
+            // pass nullptr to pdgemm_ and crash with SIGSEGV.
+            //
+            // Until the legacy path is restored or re-implemented, we
+            // explicitly reject dft_plus_u==2 with cal_force or cal_stress
+            // enabled. SCF-only runs (no force/stress) are unaffected
+            // because the energy is computed in cal_energy_correction,
+            // which does not touch DSloc arrays. Use dft_plus_u=1 for
+            // force/stress calculations.
+            if (isforce || isstress)
+            {
+                ModuleBase::WARNING_QUIT("Force_Stress_LCAO::getForceStress",
+                    "dft_plus_u==2 with cal_force or cal_stress is currently broken; "
+                    "please use dft_plus_u=1 instead. See notes in source/source_lcao/force_stress_lcao.cpp.");
+            }
             ForceStressArrays fsr_dftu;
             std::vector<std::vector<double>>* dmk_d = nullptr;
             std::vector<std::vector<std::complex<double>>>* dmk_c = nullptr;

source/source_lcao/force_stress_lcao.h

@@ -16,7 +16,7 @@
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
 #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251010
 #include "source_lcao/setup_dm.h" // mohan add 2025-11-03
-#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/exx_info.h"
 
 namespace vdw

source/source_lcao/hamilt_lcao.cpp

@@ -3,7 +3,7 @@
 #include "source_base/global_variable.h"
 #include "source_base/memory_recorder.h"
 #include "source_base/timer.h"
-#include "source_lcao/module_dftu/dftu.h"
+#include "source_lcao/module_dftu/dftu_lcao.h"
 #include "source_lcao/setup_exx.h"
 #include "source_lcao/setup_deepks.h"
 #include "source_estate/module_dm/density_matrix.h"
@@ -31,12 +31,12 @@
 #include "source_lcao/module_deltaspin/spin_constrain.h"
 #include "source_hamilt/module_hcontainer/hcontainer_funcs.h"
 #include "source_hsolver/hsolver_lcao.h"
-#include "module_dftu/dftu_lcao.h"
+#include "module_dftu/dftu_lcao_op.h"
 #include "module_operator_lcao/dspin_lcao.h"
 #include "module_operator_lcao/ekinetic.h"
 #include "module_operator_lcao/meta_lcao.h"
 #include "module_operator_lcao/nonlocal.h"
-#include "module_dftu/op_dftu_lcao.h"
+#include "module_dftu/dftu_lcao_op_legacy.h"
 #include "module_operator_lcao/op_exx_lcao.h"
 #include "module_operator_lcao/overlap.h"
 #include "module_operator_lcao/td_ekinetic_lcao.h"

source/source_lcao/hamilt_lcao.h

@@ -22,7 +22,7 @@ namespace elecstate { template <typename TK, typename TR> class DensityMatrix; }
 // Setup_DeePKS forward declaration, full definition in setup_deepks.h (moved to .cpp)
 // mohan add 20260605
 template <typename TK> class Setup_DeePKS;
-// Plus_U forward declaration, full definition in module_dftu/dftu.h (moved to .cpp)
+// Plus_U forward declaration, full definition in module_dftu/dftu_lcao.h (moved to .cpp)
 // mohan add 20260605
 class Plus_U;