Refactor DFT+U step 2

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_operator_lcao/op_dftu_lcao.h"
+#include "source_lcao/module_dftu/op_dftu_lcao.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_operator_lcao/op_dftu_lcao.h"
+#include "source_lcao/module_dftu/op_dftu_lcao.h"
 #include "source_lcao/module_operator_lcao/veff_lcao.h"
 #include "source_lcao/spar_hsr.h"
 #ifdef __EXX

source/source_lcao/CMakeLists.txt

@@ -15,7 +15,6 @@ if(ENABLE_LCAO)
         module_operator_lcao/veff_lcao.cpp
         module_operator_lcao/veff_dh.cpp
         module_operator_lcao/meta_lcao.cpp
-        module_operator_lcao/op_dftu_lcao.cpp
         module_operator_lcao/deepks_lcao.cpp
         module_operator_lcao/op_exx_lcao.cpp
         module_operator_lcao/overlap.cpp
@@ -32,8 +31,6 @@ if(ENABLE_LCAO)
         module_operator_lcao/td_pot_hybrid_fs.cpp
         module_operator_lcao/dspin_lcao.cpp
         module_operator_lcao/dspin_fs.cpp
-        module_operator_lcao/dftu_lcao.cpp
-        module_operator_lcao/dftu_fs.cpp
         module_operator_lcao/operator_fs_utils.cpp
         setup_dftu_lcao.cpp
         pulay_fs_center2.cpp

source/source_lcao/force_stress_lcao.cpp

@@ -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_operator_lcao/dftu_lcao.h"
+#include "source_lcao/module_dftu/dftu_lcao.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"

source/source_lcao/hamilt_lcao.cpp

@@ -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_operator_lcao/dftu_lcao.h"
+#include "module_dftu/dftu_lcao.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_operator_lcao/op_dftu_lcao.h"
+#include "module_dftu/op_dftu_lcao.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/module_dftu/CMakeLists.txt

@@ -8,6 +8,14 @@ list(APPEND objects
       dftu_hamilt.cpp
 )
 
+if(ENABLE_LCAO)
+  list(APPEND objects
+        dftu_lcao.cpp
+        dftu_fs.cpp
+        op_dftu_lcao.cpp
+  )
+endif()
+
 add_library(
     dftu
     OBJECT

source/source_lcao/module_dftu/dftu.cpp

@@ -105,7 +105,7 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell,
 {
     ModuleBase::TITLE("Plus_U", "cal_energy_correction");
     ModuleBase::timer::start("Plus_U", "cal_energy_correction");
-    if (!is_locale_initialized())
+    if (!is_occ_mat_initialized())
     {
         ModuleBase::timer::end("Plus_U", "cal_energy_correction");
         return;
@@ -158,11 +158,11 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell,
 
                             for (int m0 = 0; m0 < 2 * l + 1; m0++)
                             {
-                                nm_trace += this->locale[iat][l][n][spin](m0, m0);
+                                nm_trace += this->occ_mat[iat][l][n][spin](m0, m0);
                                 for (int m1 = 0; m1 < 2 * l + 1; m1++)
                                 {
-                                    nm2_trace += this->locale[iat][l][n][spin](m0, m1)
-                                                 * this->locale[iat][l][n][spin](m1, m0);
+                                    nm2_trace += this->occ_mat[iat][l][n][spin](m0, m1)
+                                                 * this->occ_mat[iat][l][n][spin](m1, m0);
                                 }
                             }
                             if (Yukawa)
@@ -186,16 +186,16 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell,
                             for (int ipol0 = 0; ipol0 < this->npol; ipol0++)
                             {
                                 const int m0_all = m0 + (2 * l + 1) * ipol0;
-                                nm_trace += this->locale[iat][l][n][0](m0_all, m0_all);
+                                nm_trace += this->occ_mat[iat][l][n][0](m0_all, m0_all);
 
                                 for (int m1 = 0; m1 < 2 * l + 1; m1++)
                                 {
                                     for (int ipol1 = 0; ipol1 < this->npol; ipol1++)
                                     {
                                         int m1_all = m1 + (2 * l + 1) * ipol1;
 
-                                        nm2_trace += this->locale[iat][l][n][0](m0_all, m1_all)
-                                                     * this->locale[iat][l][n][0](m1_all, m0_all);
+                                        nm2_trace += this->occ_mat[iat][l][n][0](m0_all, m1_all)
+                                                     * this->occ_mat[iat][l][n][0](m1_all, m0_all);
                                     }
                                 }
                             }
@@ -228,14 +228,14 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell,
                                         {
                                             double VU = 0.0;
                                             VU = get_onebody_eff_pot(T, iat, l, n, is, m1_all, m2_all, false);
-                                            energy_dc += VU * this->locale[iat][l][n][is](m1_all, m2_all);
+                                            energy_dc += VU * this->occ_mat[iat][l][n][is](m1_all, m2_all);
                                         }
                                     }
                                     else if (Plus_U::nspin == 4)
                                     {
                                         double VU = 0.0;
                                         VU = get_onebody_eff_pot(T, iat, l, n, 0, m1_all, m2_all, false);
-                                        energy_dc += VU * this->locale[iat][l][n][0](m1_all, m2_all);
+                                        energy_dc += VU * this->occ_mat[iat][l][n][0](m1_all, m2_all);
                                     }
                                 }
                             }

source/source_lcao/module_operator_lcao/dftu_fs.cpp → source/source_lcao/module_dftu/dftu_fs.cpp

@@ -138,7 +138,7 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
         }
         // first iteration to calculate occupation matrix
         std::vector<double> occ(tlp1 * tlp1 * this->nspin, 0);
-        this->dftu->get_locale_flat(iat0, target_L, occ);
+        this->dftu->get_occ_mat_flat(iat0, target_L, occ);
 
         // calculate VU
         const double u_value = this->dftu->U[T0];

source/source_lcao/module_dftu/dftu_hamilt.cpp

@@ -12,7 +12,7 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik,
 		const int npol)
 {
     ModuleBase::TITLE("Plus_U", "cal_eff_pot_c");
-    if (!is_locale_initialized())
+    if (!is_occ_mat_initialized())
     {
         return;
     }
@@ -65,7 +65,7 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik,
 void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector<int>& isk, const double* sk, const int npol)
 {
     ModuleBase::TITLE("Plus_U", "cal_eff_pot_r");
-    if (!is_locale_initialized())
+    if (!is_occ_mat_initialized())
     {
         return;
     }

source/source_lcao/module_operator_lcao/dftu_lcao.cpp → source/source_lcao/module_dftu/dftu_lcao.cpp

@@ -177,25 +177,25 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_nlm_all(const Parallel_Orbi
  * @brief Contribute DFT+U Hamiltonian to real-space HR matrix
  * 
  * @details This function handles different scenarios based on:
- * 1. Whether locale (occupation matrix) is read from file (is_locale_initialized)
+ * 1. Whether occ_mat (occupation matrix) is read from file (is_occ_mat_initialized)
  * 2. Spin configuration (nspin=1, 2, or 4)
  * 3. SCF iteration stage (first vs subsequent iterations)
  * 
- * Case 1: Locale NOT initialized (!is_locale_initialized)
+ * Case 1: Occ_mat NOT initialized (!is_occ_mat_initialized)
  *   - First electronic iteration: calculates occupation matrix from density matrix (DMR)
  *     * Uses get_dmr(current_spin) to get real-space density matrix
  *     * Accumulates contributions from all atom pairs via cal_occ()
  *     * Performs MPI reduction to sum occ across processes
- *     * Stores result via set_locale_flat() for use in VU calculation
+ *     * Stores result via set_occ_mat_flat() for use in VU calculation
  *     * For nspin=1: occ is scaled by 0.5 (since only one spin channel computed)
- *   - Subsequent iterations: locale is computed fresh each iteration from updated DMR
+ *   - Subsequent iterations: occ_mat is computed fresh each iteration from updated DMR
  * 
- * Case 2: Locale IS initialized (is_locale_initialized, i.e., read from dm_onsite.txt file)
- *   - First electronic iteration: uses pre-read locale directly without DMR calculation
+ * Case 2: Occ_mat IS initialized (is_occ_mat_initialized, i.e., read from dm_onsite.txt file)
+ *   - First electronic iteration: uses pre-read occ_mat directly without DMR calculation
  *     * Skips DMR-based occ calculation entirely
- *     * Reads locale from stored data via get_locale()
+ *     * Reads locale from stored data via get_occ_mat()
  *     * Different indexing for nspin=4 vs nspin=1/2 (see below)
- *   - After first iteration: mark_locale_dirty() is called to force recomputation
+ *   - After first iteration: mark_occ_mat_dirty() is called to force recomputation
  * 
  * Spin configurations:
  *   nspin=1 (non-spin-polarized):
@@ -206,14 +206,14 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_nlm_all(const Parallel_Orbi
  *   nspin=2 (collinear spin-polarized):
  *     - Two separate spin channels (spin-up: 0, spin-down: 1)
  *     - current_spin toggles between 0 and 1 across iterations
- *     - mark_locale_dirty() called when current_spin == 1 (last spin)
+ *     - mark_occ_mat_dirty() called when current_spin == 1 (last spin)
  *     - HR accumulated separately for each spin
  *   
  *   nspin=4 (non-collinear/SOC):
  *     - Single 4x4 Pauli matrix representation per atom
  *     - occ has 4*(2l+1)^2 elements (spin_fold=4)
- *     - get_locale uses spin=0, ipol indices for Pauli blocks
- *     - mark_locale_dirty() always called (current_spin check always true)
+ *     - get_occ_mat uses spin=0, ipol indices for Pauli blocks
+ *     - mark_occ_mat_dirty() always called (current_spin check always true)
  *     - No current_spin toggling (all spins handled simultaneously)
  * 
  * @warning THREAD SAFETY: cal_HR_IJR() updates shared HR matrix entries.
@@ -228,12 +228,12 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
     ModuleBase::TITLE("DFTU", "contributeHR");
     // Early exit conditions:
     // - get_dmr(0) == nullptr: DMR not available (typical in first iteration without file input)
-    // - !is_locale_initialized(): locale not read from file AND not yet computed from DMR
+    // - !is_occ_mat_initialized(): occ_mat not read from file AND not yet computed from DMR
     // When both true, skip DFT+U contribution entirely (first iteration, no file input)
     const bool dmr_null = (this->dftu->get_dmr(0) == nullptr);
-    const bool locale_not_init = !this->dftu->is_locale_initialized();
+    const bool occ_mat_not_init = !this->dftu->is_occ_mat_initialized();
 
-    if (dmr_null && locale_not_init)
+    if (dmr_null && occ_mat_not_init)
     {
         return;
     }
@@ -277,13 +277,13 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
         std::vector<double> occ(tlp1 * tlp1 * spin_fold, 0.0);
 
         // ============================================================
-        // BRANCH 1: Locale NOT initialized (compute from DMR)
+        // BRANCH 1: Occ_mat NOT initialized (compute from DMR)
         // ============================================================
         // This branch is taken when:
-        // - is_locale_initialized() == false (no file read or omc != 0)
+        // - is_occ_mat_initialized() == false (no file read or omc != 0)
         // - DMR is available (get_dmr() != nullptr)
         // Typical scenario: normal SCF iterations after first step
-        if (!this->dftu->is_locale_initialized())
+        if (!this->dftu->is_occ_mat_initialized())
         {
             // TODO: UNSAFE - get_dmr(current_spin) assumes DMR has correct spin indexing.
             // For nspin=2, current_spin must be correctly toggled (0 then 1).
@@ -328,37 +328,37 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
             {
                 for (auto& v : occ) { v *= 0.5; }
             }
-            this->dftu->set_locale_flat(iat0, target_L, this->current_spin, occ);
+            this->dftu->set_occ_mat_flat(iat0, target_L, this->current_spin, occ);
         }
         // ============================================================
-        // BRANCH 2: Locale IS initialized (use pre-read data)
+        // BRANCH 2: Occ_mat IS initialized (use pre-read data)
         // ============================================================
         // This branch is taken when:
-        // - is_locale_initialized() == true (locale read from dm_onsite.txt file)
+        // - is_occ_mat_initialized() == true (occ_mat read from dm_onsite.txt file)
         // - OR omc != 0 (occupation matrix control with dm_onsite_ini.txt)
         // Typical scenario: first SCF iteration with file input, or restart calculation
         else
         {
             // nspin=4: Non-collinear case with Pauli matrix representation
-            // Locale stored as single 4x4 block per atom, with spin indices embedded
+            // Occ_mat stored as single 4x4 block per atom, with spin indices embedded
             // in the matrix indices (ipol0, ipol1 for Pauli block indices)
             if (this->nspin == 4)
             {
-                // For nspin=4, locale is stored as 4 stacked tlp1^2 blocks
+                // For nspin=4, occ_mat is stored as 4 stacked tlp1^2 blocks
                 // at offsets 0, tlp1^2, 2*tlp1^2, 3*tlp1^2 for the 4 Pauli channels.
-                // Use get_locale_flat to read the stacked blocks directly
-                this->dftu->get_locale_flat(iat0, target_L, occ);
+                // Use get_occ_mat_flat to read the stacked blocks directly
+                this->dftu->get_occ_mat_flat(iat0, target_L, occ);
             }
             // nspin=1 or nspin=2: Collinear spin case
-            // Locale stored separately for each spin channel
+            // Occ_mat stored separately for each spin channel
             else
             {
                 for (int i = 0; i < static_cast<int>(occ.size()); i++)
                 {
                     // TODO: UNSAFE - current_spin must be correct for nspin=2.
-                    // If current_spin is not toggled properly, wrong spin channel's locale is read.
+                    // If current_spin is not toggled properly, wrong spin channel's occ_mat is read.
                     // This can happen if contributeHR() is called out of expected order.
-                    occ[i] = this->dftu->get_locale(iat0, target_L, 0, this->current_spin,
+                    occ[i] = this->dftu->get_occ_mat(iat0, target_L, 0, this->current_spin,
                                                       i / (2 * target_L + 1), i % (2 * target_L + 1));
                 }
             }
@@ -428,26 +428,26 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
         ModuleBase::timer::end("DFTU", "cal_vu");
     }
 
-    // 6. Post-processing: Energy correction and locale state management
+    // 6. Post-processing: Energy correction and occ_mat state management
     // For nspin=1: DFT+U energy computed for single spin channel, but should count both spins
     // set_double_energy() doubles the energy to account for degenerate spin-up/down
 	if (this->nspin == 1) 
 	{
         this->dftu->set_double_energy();
 	}
 
-    // 7. Mark locale as dirty to force recomputation in next iteration
+    // 7. Mark occ_mat as dirty to force recomputation in next iteration
     // This is called when:
     // - nspin=4: Always (all spins handled simultaneously, current_spin==0==nspin-1)
     // - nspin=2: When current_spin==1 (after spin-down calculation, last spin channel)
     // - nspin=1: When current_spin==0==nspin-1 (always called)
     // 
-    // Purpose: Ensure locale is recomputed from updated DMR in next SCF iteration,
+    // Purpose: Ensure occ_mat is recomputed from updated DMR in next SCF iteration,
     // rather than using stale pre-read data from file.
     // TODO: This logic is confusing. Consider explicit variable like `is_last_spin_channel`.
 	if (this->current_spin == this->nspin - 1 || this->nspin == 4) 
 	{
-		this->dftu->mark_locale_dirty();
+		this->dftu->mark_occ_mat_dirty();
 	}
 
     // 8. Spin channel toggling for nspin=2