DFT+U refactor, step 8

source/CMakeLists.txt

@@ -695,7 +695,6 @@ target_link_libraries(
     psi_overall_init
     psi_init
     psi
-    dftu
     deltaspin
     container
     device
@@ -711,6 +710,7 @@ if(ENABLE_LCAO)
     PRIVATE
       hamilt_lcao
       tddft
+      dftu
       orb
       gint
       hcontainer

source/Makefile.Objects

@@ -837,8 +837,8 @@ OBJS_SRCPW=h_ewald_pw.o\
     update_cell_pw.o\
     dftu_base.o\
     dftu_base_io.o\
-    dftu_base_occ.o\
-    dftu_base_tools.o\
+    dftu_pw.o\
+    dftu_pw_tools.o\
     yukawa_screening.o\
     setup_dftu_pw.o\
     deltaspin_pw.o\

source/source_esolver/esolver_dfpt_pw.cpp

@@ -354,7 +354,7 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
 
     dfpt_->init(ucell, *this->stp.psi_cpu, this->pw_rho, this->pw_wfc, &this->sf, veff_r,
                 this->pelec->wg, this->pelec->ekb, xc_adapter_, nelec_, ecutwfc_,
-                dft_plus_u_ ? &this->dftu : nullptr);
+                dft_plus_u_ ? this->dftu_ : nullptr);
 }
 
 void ESolver_DFPT_PW::run_post_process(UnitCell& ucell)

source/source_esolver/esolver_double_xc.cpp

@@ -148,7 +148,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                                                              this->two_center_bundle_,
                                                              this->orb_,
                                                              this->dmat_base.dm,
-                                                             &this->dftu,
+                                                             this->dftu_,
                                                              this->deepks,
                                                              istep,
                                                              this->exx_nao,
@@ -413,7 +413,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
                        this->kv,
                        this->pw_rho,
                        this->solvent,
-                       this->dftu,
+                       *this->dftu_,
                        this->deepks,
                        this->exx_nao,
                        &ucell.symm,

source/source_esolver/esolver_fp.cpp

@@ -153,7 +153,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
     CE.update_delta_rho(ucell, &(this->chr), &(this->sf));
 
     //! print out charge density, potential, elf, etc.
-	ModuleIO::ctrl_output_fp(ucell, this->pelec, this->pw_big, this->pw_rhod, 
+	ModuleIO::ctrl_output_fp(ucell, *this->inp_, this->pelec, this->pw_big, this->pw_rhod, 
 			this->chr, this->solvent, this->Pgrid, istep); 
 
 }

source/source_esolver/esolver_ks.cpp

@@ -246,7 +246,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &
 #ifdef __LCAO
     if (this->inp_->dft_plus_u)
     {
-        converged_u = this->dftu.u_converged();
+        converged_u = this->dftu_->u_converged();
     }
 #endif
 

source/source_esolver/esolver_ks.h

@@ -60,8 +60,9 @@ class ESolver_KS : public ESolver_FP
     //! nonlocal pseudopotentials
     pseudopot_cell_vnl ppcell;
 
-    //! DFT+U method, mohan add 2025-11-07
-    Plus_U dftu;
+    //! DFT+U method (base-class pointer; PW news Plus_U_Base, LCAO news Plus_U),
+    //! allocated in the derived-class constructor, deleted in ~ESolver_KS.
+    Plus_U_Base* dftu_ = nullptr;
 
     std::string basisname;          //! esolver_ks_lcao.cpp
     double esolver_KS_ne = 0.0;     //! number of electrons

source/source_esolver/esolver_ks_lcao.cpp

@@ -5,6 +5,7 @@
 #include "source_lcao/module_deltaspin/spin_constrain.h"
 #include "source_lcao/module_deltaspin/deltaspin_lcao.h"
 #include "source_lcao/setup_dftu_lcao.h"
+#include "source_lcao/module_dftu/dftu_nao.h" // Plus_U (LCAO DFT+U derived class)
 #include "source_hamilt/hs_matrix_k.h"
 #include "source_estate/module_charge/symm_rho.h"
 #include "source_lcao/lcao_domain.h" // need DeePKS_init
@@ -35,6 +36,7 @@ ESolver_KS_LCAO<TK, TR>::ESolver_KS_LCAO()
 {
     this->classname = "ESolver_KS_LCAO";
     this->basisname = "LCAO";
+    this->dftu_ = new Plus_U();
 }
 
 template <typename TK, typename TR>
@@ -92,7 +94,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(BaseCell& basecell, const Input
       this->dmat, this->chr, inp);
 
     LCAO_domain::set_pot<TK>(ucell, this->kv, this->sf, *this->pw_rho, *this->pw_rhod,
-      this->pelec, this->orb_, this->pv, this->locpp, this->dftu,
+      this->pelec, this->orb_, this->pv, this->locpp, *this->dftu_,
       this->solvent, this->exx_nao, this->deepks, inp, this->exx_info_);
 
     //! if kpar is not divisible by nks, print a warning
@@ -156,7 +158,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     {
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
-            two_center_bundle_, orb_, this->dmat.dm, &this->dftu, this->deepks, istep, exx_nao, this->exx_info_);
+            two_center_bundle_, orb_, this->dmat.dm, this->dftu_, this->deepks, istep, exx_nao, this->exx_info_);
     }
 
     // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
@@ -256,7 +258,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
                        this->gd, this->pv, this->pelec, this->dmat, this->psi,
                        two_center_bundle_, orb_, force, this->scs,
                        this->locpp, this->sf, this->kv,
-                       this->pw_rho, this->solvent, this->dftu, this->deepks,
+                       this->pw_rho, this->solvent, *this->dftu_, this->deepks,
                        this->exx_nao, &ucell.symm, this->exx_info_, this->inp_->td_stype,
                        static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt));
 
@@ -333,7 +335,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     // call iter_init() of ESolver_KS
     ESolver_KS::iter_init(ucell, istep, iter);
 
-    module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dftu, 
+    module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, *this->dftu_, 
       this->dmat.dm->get_DMR_pointer(1)->get_nnr(), *this->inp_); 
 
     if (iter == 1)
@@ -387,7 +389,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     }
 #endif
 
-    init_dftu_lcao<TK>(istep, iter, this->inp_->dft_plus_u, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx);
+    init_dftu_lcao<TK>(istep, iter, this->inp_->dft_plus_u, this->dftu_, this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx);
 
 #ifdef __MLALGO
     // the density matrixes of DeePKS have been updated in each iter
@@ -508,14 +510,14 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 	const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_DMK_vector();
 
     // 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, PARAM.globalv.gamma_only_local);
+    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, PARAM.globalv.gamma_only_local);
 
     // 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_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());
+        this->pelec->set_dftu_energy(this->dftu_->get_energy());
     }
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
@@ -579,7 +581,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
     //! 2) output of lcao every few ionic steps
     ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
             *this->inp_, this->kv, this->pelec, this->dmat.dm, this->pv,
-            this->gd, this->psi, hamilt_lcao, this->dftu, this->two_center_bundle_,
+            this->gd, this->psi, hamilt_lcao, *this->dftu_, this->two_center_bundle_,
             this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
             this->pw_rhod, this->locpp.vloc, this->solvent,
             this->rdmft_solver, this->deepks, this->exx_nao, this->exx_info_,

source/source_esolver/esolver_ks_pw.cpp

@@ -32,6 +32,8 @@ ESolver_KS_PW<T, Device>::ESolver_KS_PW()
 {
     this->classname = "ESolver_KS_PW";
     this->basisname = "PW";
+    // PW basis: the DFT+U object is the base class (no LCAO orbitals).
+    this->dftu_ = new Plus_U_Base();
 }
 
 template <typename T, typename Device>
@@ -65,7 +67,7 @@ void ESolver_KS_PW<T, Device>::allocate_hamilt(const UnitCell& ucell)
                                                      this->pw_wfc,
                                                      &this->kv,
                                                      &this->ppcell,
-                                                     &this->dftu,
+                                                     this->dftu_,
                                                      &ucell,
                                                      &this->general_exx_info_);
 }
@@ -171,7 +173,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
                   this->chr,
                   this->locpp,
                   this->ppcell,
-                  this->dftu,
+                  *this->dftu_,
                   this->vsep_cell,
                   this->stp.template get_psi_t<T, Device>(),
                   this->p_hamilt,
@@ -194,7 +196,7 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
 {
     ESolver_KS::iter_init(ucell, istep, iter);
 
-    module_charge::chgmixing_ks_pw(iter, this->p_chgmix, this->dftu, *this->inp_);
+    module_charge::chgmixing_ks_pw(iter, this->p_chgmix, *this->dftu_, *this->inp_);
 
     // mohan move harris functional here, 2012-06-05
     // use 'rho(in)' and 'v_h and v_xc'(in)
@@ -204,7 +206,7 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
     // should before lambda loop in DeltaSpin
     DFTU_BASE::iter_init_dftu_pw(iter,
                           istep,
-                          this->dftu,
+                          *this->dftu_,
                           this->stp.template get_psi_t<T, Device>(),
                           this->pelec->wg,
                           ucell,
@@ -214,7 +216,7 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
     // 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());
+        this->pelec->set_dftu_energy(this->dftu_->get_energy());
     }
 }
 
@@ -383,7 +385,7 @@ void ESolver_KS_PW<T, Device>::cal_force(BaseCell& basecell, ModuleBase::matrix&
                  &ucell.symm,
                  &this->sf,
                  this->solvent,
-                 &this->dftu,
+                 this->dftu_,
                  &this->locpp,
                  &this->ppcell,
                  &this->kv,
@@ -405,7 +407,7 @@ void ESolver_KS_PW<T, Device>::cal_stress(BaseCell& basecell, ModuleBase::matrix
     ss.cal_stress(stress,
                   ucell,
                   this->get_vdw_result(),
-                  this->dftu,
+                  *this->dftu_,
                   this->locpp,
                   this->ppcell,
                   this->pw_rhod,

source/source_esolver/lcao_others.cpp

@@ -133,7 +133,7 @@ void ESolver_KS_LCAO<TK, TR>::others(BaseCell& basecell, const int istep)
                                                         two_center_bundle_,
                                                         orb_,
                                                         this->dmat.dm,
-                                                        &this->dftu,
+                                                        this->dftu_,
                                                         this->deepks,
                                                         istep,
                                                         this->exx_nao,

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
-    // (pot_uterm_pw.size() = pot_index * 2 for nspin=2)
+    // (uterm_mat.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/occ_matrix.cpp

@@ -4,7 +4,7 @@
 #include "source_cell/unitcell.h"
 
 void OccupationMatrix::init(const UnitCell& cell,
-                            const std::vector<int>& orbital_corr,
+                            const std::vector<int>& l_channel,
                             const int nspin,
                             const int npol)
 {
@@ -25,7 +25,7 @@ void OccupationMatrix::init(const UnitCell& cell,
             occ_save_[iat].resize(cell.atoms[it].nwl + 1);
             iatlnmipol2iwt_[iat].resize(cell.atoms[it].nwl + 1);
 
-            if (orbital_corr[it] == -1)
+            if (l_channel[it] == -1)
             {
                 continue;
             }
@@ -123,11 +123,11 @@ void OccupationMatrix::set_flat(const int iat, const int l, const int spin,
     }
 }
 
-void OccupationMatrix::zero(const UnitCell& cell, const std::vector<int>& orbital_corr)
+void OccupationMatrix::zero(const UnitCell& cell, const std::vector<int>& l_channel)
 {
     for (int T = 0; T < cell.ntype; T++)
     {
-        if (orbital_corr[T] == -1)
+        if (l_channel[T] == -1)
         {
             continue;
         }
@@ -157,14 +157,14 @@ void OccupationMatrix::zero(const UnitCell& cell, const std::vector<int>& orbita
     }
 }
 
-void OccupationMatrix::copy_to_save(const UnitCell& cell, const std::vector<int>& orbital_corr)
+void OccupationMatrix::copy_to_save(const UnitCell& cell, const std::vector<int>& l_channel)
 {
     ModuleBase::TITLE("OccupationMatrix", "copy_to_save");
     ModuleBase::timer::start("OccupationMatrix", "copy_to_save");
 
     for (int T = 0; T < cell.ntype; T++)
     {
-        const int target_l = orbital_corr[T];
+        const int target_l = l_channel[T];
         if (target_l == -1)
         {
             continue;
@@ -189,7 +189,7 @@ void OccupationMatrix::copy_to_save(const UnitCell& cell, const std::vector<int>
 }
 
 void OccupationMatrix::write_to_flat(const UnitCell& cell,
-                                     const std::vector<int>& orbital_corr,
+                                     const std::vector<int>& l_channel,
                                      const std::vector<int>& index,
                                      std::vector<double>& uom) const
 {
@@ -200,7 +200,7 @@ void OccupationMatrix::write_to_flat(const UnitCell& cell,
     for (int iat = 0; iat < cell.nat; iat++)
     {
         const int it = cell.iat2it[iat];
-        const int target_l = orbital_corr[it];
+        const int target_l = l_channel[it];
         if (target_l == -1)
         {
             continue;
@@ -223,13 +223,13 @@ void OccupationMatrix::write_to_flat(const UnitCell& cell,
 }
 
 void OccupationMatrix::read_from_flat(const UnitCell& cell,
-                                      const std::vector<int>& orbital_corr,
+                                      const std::vector<int>& l_channel,
                                       const std::vector<int>& index,
                                       const std::vector<double>& uom)
 {
     for (int T = 0; T < cell.ntype; T++)
     {
-        const int l = orbital_corr[T];
+        const int l = l_channel[T];
         if (l == -1)
         {
             continue;
@@ -261,7 +261,7 @@ void OccupationMatrix::read_from_flat(const UnitCell& cell,
 }
 
 void OccupationMatrix::write_save_to_flat(const UnitCell& cell,
-                                          const std::vector<int>& orbital_corr,
+                                          const std::vector<int>& l_channel,
                                           const std::vector<int>& index,
                                           std::vector<double>& uom_save) const
 {
@@ -271,7 +271,7 @@ void OccupationMatrix::write_save_to_flat(const UnitCell& cell,
     }
     for (int T = 0; T < cell.ntype; T++)
     {
-        const int target_l = orbital_corr[T];
+        const int target_l = l_channel[T];
         if (target_l == -1)
         {
             continue;
@@ -318,13 +318,13 @@ namespace elecstate
 void mix_occ_with_save(std::vector<std::vector<std::vector<std::vector<ModuleBase::matrix>>>>& occ_mat,
                        const std::vector<std::vector<std::vector<std::vector<ModuleBase::matrix>>>>& occ_mat_save,
                        const UnitCell& cell,
-                       const std::vector<int>& orbital_corr,
+                       const std::vector<int>& l_channel,
                        const int nspin,
                        const double beta)
 {
     for (int T = 0; T < cell.ntype; T++)
     {
-        const int target_l = orbital_corr[T];
+        const int target_l = l_channel[T];
         if (target_l == -1)
         {
             continue;