Refactor DFT+U, step 4

source/Makefile.Objects

@@ -742,7 +742,8 @@ OBJS_SRCPW=h_ewald_pw.o\
     update_cell_pw.o\
     dftu_base.o\
     dftu_output.o\
-    dftu_pw.o\
+    dftu_tools_pw.o\
+    dftu_cal_occ_pw.o\
     setup_dftu_pw.o\
     deltaspin_pw.o\
     force_pw.o\

source/source_esolver/esolver_ks_lcao.cpp

@@ -510,7 +510,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 
     // 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).
+    // via cal_pot_onsite) 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());

source/source_estate/module_charge/charge_mixing.cpp

@@ -264,7 +264,7 @@ void Charge_Mixing::allocate_mixing_uom(int uom_size)
     ModuleBase::TITLE("Charge_Mixing", "allocate_mixing_uom");
     ModuleBase::timer::start("Charge_Mixing", "allocate_mixing_uom");
     // For nspin=2, uom_size already includes both spin channels
-    // (eff_pot_pw.size() = pot_index * 2 for nspin=2)
+    // (pot_uterm_pw.size() = pot_index * 2 for nspin=2)
     // So uom_fold should always be 1
     this->mixing->init_mixing_data(this->uom_mdata, uom_size, sizeof(double));
     this->uom_mdata.reset();

source/source_estate/module_charge/chgmixing.cpp

@@ -133,7 +133,7 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
             // enable mixing_dftu for DFT+U occupation mixing
             dftu.enable_mixing();
             // allocate memory for uom_mdata
-            p_chgmix->allocate_mixing_uom(dftu.get_size_eff_pot_pw());
+            p_chgmix->allocate_mixing_uom(dftu.get_size_pot_uterm_pw());
         }
     }
 
@@ -145,11 +145,11 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
 
         if (inp.dft_plus_u)
         {
-            if (dftu.uramping > 0.01 && !dftu.u_converged())
+            if (dftu.get_uramping() > 0.01 && !dftu.u_converged())
             {
                 p_chgmix->mixing_restart_step = inp.scf_nmax + 1;
             }
-            if (dftu.uramping > 0.01)
+            if (dftu.get_uramping() > 0.01)
             {
                 bool do_uramping = true;
                 if (inp.sc_mag_switch)
@@ -197,7 +197,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number
             dftu.enable_mixing();
         }
         // this output will be removed once the feeature is stable
-        if (dftu.uramping > 0.01)
+        if (dftu.get_uramping() > 0.01)
         {
             std::cout << " U-Ramping! Current U = ";
             for (int i = 0; i < dftu.get_num_u_types(); i++)
@@ -216,7 +216,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number
         if (inp.dft_plus_u)
         {
             dftu.uramping_update(); // update U by uramping if uramping > 0.01
-            if (dftu.uramping > 0.01)
+            if (dftu.get_uramping() > 0.01)
             {
                 std::cout << " U-Ramping! Current U = ";
                 for (int i = 0; i < dftu.get_num_u_types(); i++)
@@ -225,7 +225,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number
                 }
                 std::cout << " eV " << std::endl;
             }
-            if (dftu.uramping > 0.01 && !dftu.u_converged())
+            if (dftu.get_uramping() > 0.01 && !dftu.u_converged())
             {
                 p_chgmix->mixing_restart_step = inp.scf_nmax + 1;
             }

source/source_io/module_parameter/input_parameter.h

@@ -112,7 +112,7 @@ struct Input_para
     double mixing_gg0_min = 0.1;
     double mixing_angle = -10.0;
     bool mixing_tau = false;  ///< whether to mix tau in mgga
-    bool mixing_dftu = false; ///< whether to mix locale in DFT+U
+    bool mixing_dftu = false; ///< whether to mix occ_mat in DFT+U
     bool mixing_dmr = false;  ///< whether to mix real space density matrix
 
     bool gamma_only = false;   ///< for plane wave.

source/source_io/module_parameter/read_inp_estruc.cpp

@@ -757,7 +757,7 @@ This setting takes effect only when the selected exchange-correlation functional
     }
     {
         Input_Item item("mixing_dftu");
-        item.annotation = "whether to mix locale in DFT+U calculation";
+        item.annotation = "whether to mix occ_mat in DFT+U calculation";
         item.category = "Electronic structure";
         item.type = "Boolean";
         item.description = R"(Whether to mix the occupation matrices.

source/source_lcao/force_stress_lcao.cpp

@@ -2,6 +2,7 @@
 
 #include "source_base/parallel_reduce.h"
 #include "source_lcao/module_dftu/dftu_lcao.h" //Quxin add for DFT+U on 20201029
+#include "source_lcao/module_dftu/dftu_force.h"
 #include "source_io/module_output/output_log.h"
 #include "source_io/module_parameter/parameter.h"
 // new
@@ -456,7 +457,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
             std::vector<std::vector<double>>* dmk_d = nullptr;
             std::vector<std::vector<std::complex<double>>>* dmk_c = nullptr;
             assign_dmk_ptr<T>(dmat.dm, dmk_d, dmk_c, PARAM.globalv.gamma_only_local);
-            dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv, PARAM.globalv.npol);
+            DFTU_LCAO::force_stress(dftu, ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv, PARAM.globalv.npol);
         }
         else
         {