Merge branch 'develop' into 2026-08-28-2

Author: mohanchen

AGENTS.md

@@ -98,6 +98,12 @@ rules. Read the complete governance document before making or reviewing changes:
   ask the developer whether to write them in; be cautious and skip unclear
   or unverified lessons.
 
+## Refactoring Patterns
+
+- Member -> free function: inventory `this->` reads; pass as params (const
+  for config, ref for mutable state); move only when body is `this`-free;
+  keep thin wrapper; compile each step.
+
 ## Local Commands
 
 ```bash

source/Makefile.Objects

@@ -818,7 +818,7 @@ OBJS_DFTU=dftu.o\
       dftu_yukawa.o\
       dftu_folding.o\
       dftu_tools.o\
-      dftu_occup.o\
+      dftu_lcao_occ.o\
       dftu_hamilt.o\
       setup_dftu_pw.o
 

source/source_esolver/esolver_ks_lcao.cpp

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

source/source_io/module_hs/write_vxc.hpp

@@ -209,7 +209,7 @@ void write_Vxc(const int nspin,
         &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, &exx_info, hamilt::Add_Hexx_Type::k);
     std::vector<std::vector<double>> e_orb_exx; // orbital energy (EXX)
 #endif
-    hamilt::OperatorDFTU<hamilt::OperatorLCAO<TK, TR>> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, nullptr, kv.isk, PARAM.globalv.npol);
+    hamilt::OperatorDFTU<hamilt::OperatorLCAO<TK, TR>> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, ucell, nullptr, kv.isk, PARAM.globalv.npol);
 
     // 4. calculate and write the MO-matrix Exc
     Parallel_2D p2d;

source/source_lcao/force_stress_lcao.cpp

@@ -457,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_LCAO::force_stress(dftu, ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv, PARAM.globalv.npol);
+            DFTU_LCAO::force_stress(dftu, isforce, isstress, ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv, PARAM.globalv.npol, PARAM.globalv.gamma_only_local);
         }
         else
         {

source/source_lcao/hamilt_lcao.cpp

@@ -226,10 +226,11 @@ HamiltLCAO<TK, TR>::HamiltLCAO(const UnitCell& ucell,
             {
                 plus_u = new OperatorDFTU<OperatorLCAO<TK, TR>>(this->hsk,
                                                               this->kv->kvec_d,
-															  this->hR,
-															  p_dftu,
-															  this->kv->isk,
-															  PARAM.globalv.npol);
+																  this->hR,
+																  ucell,
+																  p_dftu,
+                                                              this->kv->isk,
+                                                              PARAM.globalv.npol);
             }
             else
             {
@@ -383,8 +384,9 @@ HamiltLCAO<TK, TR>::HamiltLCAO(const UnitCell& ucell,
             {
                 plus_u = new OperatorDFTU<OperatorLCAO<TK, TR>>(this->hsk,
                                                               this->kv->kvec_d,
-															  this->hR,
-															  p_dftu,
+																  this->hR,
+																  ucell,
+																  p_dftu,
                                                               this->kv->isk,
                                                               PARAM.globalv.npol);
             }

source/source_lcao/lcao_set.cpp

@@ -89,10 +89,7 @@ void LCAO_domain::set_pot(
                   PARAM.globalv.global_out_dir,
                   inp.init_chg,
                   pv.get_global_row_size(),
-                  PARAM.globalv.gamma_only_local,
                   inp.ks_solver,
-                  inp.cal_force,
-                  inp.cal_stress,
                   inp.device,
                   inp.kpar,
                   PARAM.globalv.hubbard_u,

source/source_lcao/module_dftu/CMakeLists.txt

@@ -3,8 +3,9 @@ list(APPEND objects
       dftu_force.cpp
       dftu_yukawa.cpp
       dftu_folding.cpp
-      dftu_tools.cpp
-      dftu_occup.cpp
+      dftu_lcao_pots.cpp
+      dftu_lcao_occ.cpp
+      dftu_lcao_energy.cpp
       dftu_hamilt.cpp
 )
 

source/source_lcao/module_dftu/dftu_force.cpp

@@ -2,9 +2,11 @@
 #include "dftu_force.h"
 #include "dftu_folding.h"
 #include "dftu_lcao.h"
+#include "dftu_lcao_pots.h"
 #include "source_base/global_function.h"
 #include "source_base/module_external/scalapack_connector.h"
 #include "source_base/parallel_reduce.h"
+#include "source_io/module_parameter/parameter.h"
 #include "source_base/timer.h"
 
 #include <complex>
@@ -14,6 +16,8 @@
 namespace DFTU_LCAO {
 
 void force_stress(Plus_U& dftu,
+                  const bool cal_force,
+                  const bool cal_stress,
                   const UnitCell& ucell,
                   const Grid_Driver& gd,
                   std::vector<std::vector<double>>* dmk_d,
@@ -23,7 +27,8 @@ void force_stress(Plus_U& dftu,
                   ModuleBase::matrix& force_dftu,
                   ModuleBase::matrix& stress_dftu,
                   const K_Vectors& kv,
-                  const int npol)
+                  const int npol,
+                  const bool gamma_only_local)
 {
     ModuleBase::TITLE("DFTU_LCAO", "force_stress");
     ModuleBase::timer::start("DFTU_LCAO", "force_stress");
@@ -35,16 +40,16 @@ void force_stress(Plus_U& dftu,
     // fsr_dftu is created without allocation), we fail early with a clear
     // message instead of letting pdgemm_ dereference nullptr and crash.
     // See force_stress_lcao.cpp for the historical background.
-    if (dftu.is_gamma_only_local())
+    if (gamma_only_local)
     {
-        if (dftu.is_cal_force()
+        if (cal_force
             && (fsr.DSloc_x == nullptr || fsr.DSloc_y == nullptr || fsr.DSloc_z == nullptr))
         {
             ModuleBase::WARNING_QUIT("DFTU_LCAO::force_stress",
                 "fsr.DSloc_x/y/z are nullptr in gamma_only path; the caller must allocate and fill them. "
                 "See notes in source/source_lcao/force_stress_lcao.cpp.");
         }
-        if (dftu.is_cal_stress()
+        if (cal_stress
             && (fsr.DSloc_x == nullptr || fsr.DSloc_y == nullptr || fsr.DSloc_z == nullptr
                 || fsr.DH_r == nullptr))
         {
@@ -56,14 +61,14 @@ void force_stress(Plus_U& dftu,
     }
     else
     {
-        if (dftu.is_cal_force()
+        if (cal_force
             && (fsr.DSloc_Rx == nullptr || fsr.DSloc_Ry == nullptr || fsr.DSloc_Rz == nullptr))
         {
             ModuleBase::WARNING_QUIT("DFTU_LCAO::force_stress",
                 "fsr.DSloc_Rx/Ry/Rz are nullptr in multik path; the caller must allocate and fill them. "
                 "See notes in source/source_lcao/force_stress_lcao.cpp.");
         }
-        if (dftu.is_cal_stress()
+        if (cal_stress
             && (fsr.DSloc_Rx == nullptr || fsr.DSloc_Ry == nullptr || fsr.DSloc_Rz == nullptr
                 || fsr.DH_r == nullptr))
         {
@@ -79,26 +84,26 @@ void force_stress(Plus_U& dftu,
     // explicit ic * pv.nrow + ir indices. All ks_solvers accepted by INPUT
     // validation are column-major today; abort loudly instead of silently
     // producing wrong forces/stresses if that assumption ever changes.
-    if ((dftu.is_cal_force() || dftu.is_cal_stress())
-        && !ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(dftu.get_ks_solver()))
+    if ((cal_force || cal_stress)
+        && !ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
     {
         ModuleBase::WARNING_QUIT("DFTU_LCAO::force_stress",
             "non column-major ks_solver is not supported for DFT+U force/stress; "
             "the folded matrix layout assumption would be violated");
     }
 
-    const int nlocal = dftu.get_nlocal();
+    const int nlocal = pv.get_global_row_size();
 
-    if (dftu.is_cal_force())
+    if (cal_force)
     {
         force_dftu.zero_out();
     }
-    if (dftu.is_cal_stress())
+    if (cal_stress)
     {
         stress_dftu.zero_out();
     }
 
-    if (dftu.is_gamma_only_local())
+    if (gamma_only_local)
     {
         const char transN = 'N';
         const char transT = 'T';
@@ -115,7 +120,7 @@ void force_stress(Plus_U& dftu,
 
             double* pot_onsite = new double[pv.nloc];
 
-            dftu.pot_onsite_real(spin, false, pot_onsite, npol);
+            DFTU_LCAO::pot_onsite_real(dftu, ucell, &pv, spin, false, pot_onsite, npol);
 
 #ifdef __MPI
             ScalapackConnector::gemm(transT, transN, nlocal, nlocal, nlocal,
@@ -127,18 +132,18 @@ void force_stress(Plus_U& dftu,
 
             delete[] pot_onsite;
 
-            if (dftu.is_cal_force())
+            if (cal_force)
             {
-                cal_force_gamma(dftu.get_nlocal(), dftu.get_npol(),
+                cal_force_gamma(nlocal, npol,
                                 dftu.get_orbital_corr_vec(), dftu.get_iatlnmipol2iwt(),
                                 ucell, &rho_pot_onsite[0], pv,
                                 fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z, force_dftu);
             }
 
-            if (dftu.is_cal_stress())
+            if (cal_stress)
             {
-                cal_stress_gamma(dftu.get_nlocal(), dftu.get_npol(),
-                                 dftu.get_ks_solver(), dftu.get_orb_cutoff(),
+                cal_stress_gamma(nlocal, npol,
+                                 PARAM.inp.ks_solver, dftu.get_orb_cutoff(),
                                  ucell, pv, &gd,
                                  fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z, fsr.DH_r,
                                  &rho_pot_onsite[0], stress_dftu);
@@ -161,7 +166,7 @@ void force_stress(Plus_U& dftu,
 
             std::complex<double>* pot_onsite = new std::complex<double>[pv.nloc];
 
-            dftu.pot_onsite_complex(spin, false, pot_onsite, npol);
+            DFTU_LCAO::pot_onsite_complex(dftu, ucell, &pv, spin, false, pot_onsite, npol);
 
 
 #ifdef __MPI
@@ -173,28 +178,28 @@ void force_stress(Plus_U& dftu,
 
             delete[] pot_onsite;
 
-            if (dftu.is_cal_force())
+            if (cal_force)
             {
-                cal_force_k(dftu.get_nlocal(), dftu.get_npol(),
-                            dftu.get_ks_solver(), dftu.get_orb_cutoff(),
+                cal_force_k(nlocal, npol,
+                            PARAM.inp.ks_solver, dftu.get_orb_cutoff(),
                             dftu.get_orbital_corr_vec(), dftu.get_iatlnmipol2iwt(),
                             ucell, gd, fsr, pv, ik, &rho_pot_onsite[0], force_dftu, kv.kvec_d[ik]);
             }
-            if (dftu.is_cal_stress())
+            if (cal_stress)
             {
-                cal_stress_k(dftu.get_nlocal(), dftu.get_npol(),
-                             dftu.get_ks_solver(), dftu.get_orb_cutoff(),
+                cal_stress_k(nlocal, npol,
+                             PARAM.inp.ks_solver, dftu.get_orb_cutoff(),
                              ucell, gd, fsr, pv, ik, &rho_pot_onsite[0], stress_dftu, kv.kvec_d[ik]);
             }
         } // ik
     }
 
-    if (dftu.is_cal_force())
+    if (cal_force)
     {
         Parallel_Reduce::reduce_pool(force_dftu.c, force_dftu.nr * force_dftu.nc);
     }
 
-    if (dftu.is_cal_stress())
+    if (cal_stress)
     {
         Parallel_Reduce::reduce_pool(stress_dftu.c, stress_dftu.nr * stress_dftu.nc);
 

source/source_lcao/module_dftu/dftu_force.h

@@ -30,6 +30,8 @@ namespace DFTU_LCAO {
 /// Takes Plus_U& because it calls dftu.pot_onsite_real/complex,
 /// which are still members of Plus_U (defined in dftu_tools.cpp).
 void force_stress(Plus_U& dftu,
+                  const bool cal_force,
+                  const bool cal_stress,
                   const UnitCell& ucell,
                   const Grid_Driver& gd,
                   std::vector<std::vector<double>>* dmk_d,
@@ -39,7 +41,8 @@ void force_stress(Plus_U& dftu,
                   ModuleBase::matrix& force_dftu,
                   ModuleBase::matrix& stress_dftu,
                   const K_Vectors& kv,
-                  const int npol);
+                  const int npol,
+                  const bool gamma_only_local);
 
 /// @brief Force contribution at a k-point (multik path).
 void cal_force_k(int nlocal,