Merge branch 'develop' into 2026-09-02-line2-DFPT

Author: mohanchen

source/Makefile.Objects

@@ -827,7 +827,6 @@ OBJS_SRCPW=h_ewald_pw.o\
     yukawa_screening.o\
     setup_dftu_pw.o\
     deltaspin_pw.o\
-    deltaspin_pw_impl.o\
     force_pw.o\
     force_pw_us.o\
     force_pw_nl.o\
@@ -914,12 +913,14 @@ OBJS_DFTU=dftu_nao.o\
 OBJS_DELTASPIN=basic_funcs.o\
       cal_mw_from_lambda.o\
       cal_mw.o\
-      init_sc.o\
+      deltaspin_init.o\
       lambda_loop_helper.o\
       lambda_loop.o\
       spin_constrain.o\
-      cal_mw_helper.o\
       deltaspin_lcao.o\
+      deltaspin_lcao_mi.o\
+      deltaspin_state.o\
+      deltaspin_pw_mi.o\
       mi_tools.o\
       template_helpers.o\
 

source/source_estate/test/CMakeLists.txt

@@ -54,6 +54,7 @@ AddTest(
     ../occupy.cpp
     ../module_charge/charge_mpi.cpp
     ../../source_lcao/module_deltaspin/spin_constrain.cpp
+    ../../source_lcao/module_deltaspin/deltaspin_state.cpp
     ../../source_psi/psi.cpp
     ../../source_base/module_device/memory_op.cpp
     ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp

source/source_io/module_ctrl/ctrl_output_pw.cpp

@@ -11,6 +11,7 @@
 #include "../module_wf/write_wfc_pw.h" // use write_wfc_pw
 #include "source_base/formatter.h"
 #include "source_lcao/module_deltaspin/lambda_loop_helper.h"
+#include "source_lcao/module_deltaspin/deltaspin_pw_mi.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
 #include "source_pw/module_pwdft/elecond.h"
 #include "source_pw/module_pwdft/onsite_proj.h" // use projector
@@ -219,7 +220,7 @@ void ModuleIO::ctrl_scf_pw(const int istep,
     if (inp.sc_mag_switch)
     {
         spinconstrain::SpinConstrain<std::complex<double>>& sc = spinconstrain::SpinConstrain<std::complex<double>>::getScInstance();
-        sc.cal_mi_pw();
+        spinconstrain::pw::cal_mi_pw(sc.state_, sc.psi, sc.pelec);
         spinconstrain::print_Mag_Force(sc, GlobalV::ofs_running);
     }
 

source/source_lcao/module_deltaspin/CMakeLists.txt

@@ -1,14 +1,16 @@
 list(APPEND objects
     spin_constrain.cpp
-    init_sc.cpp
+    deltaspin_init.cpp
     cal_mw.cpp
     basic_funcs.cpp
     lambda_loop_helper.cpp
     lambda_loop.cpp
     cal_mw_from_lambda.cpp
     template_helpers.cpp
     deltaspin_lcao.cpp
-    cal_mw_helper.cpp
+    deltaspin_lcao_mi.cpp
+    deltaspin_state.cpp
+    deltaspin_pw_mi.cpp
     mi_tools.cpp
 )
 

source/source_lcao/module_deltaspin/cal_mw.cpp

@@ -1,103 +1,40 @@
-#include "source_base/tool_title.h"
-#include "source_base/timer.h"
-#include "spin_constrain.h"
 #ifdef __LCAO
-#include "source_estate/elecstate_lcao.h"
-#include "source_lcao/hamilt_lcao.h"
-#include "source_lcao/module_operator_lcao/dspin_lcao.h"
 
 /**
  * @file cal_mw.cpp
- * @brief Magnetic moment calculation for LCAO and PW basis sets.
- *
- * @par cal_mi_lcao (LCAO)
- * Uses the DeltaSpin operator to compute magnetic moments from the density
- * matrix via real-space projection. For nspin=2, only the z-component is
- * extracted. For nspin=4, all three components are extracted from the
- * interleaved 4-component spinor density matrix.
+ * @brief Thin LCAO shells on SpinConstrain: cal_mi_lcao() and set_operator().
  *
- * @par cal_mi_pw (PW)
- * Uses the OnsiteProjector to compute atomic projections <alpha_{l,m}|psi_{k,i}>
- * (becp coefficients), then decomposes these into magnetic moments using
- * Pauli matrix traces (accumulate_Mi_from_becp).
- *
- * @par Error conditions
- * - Dynamic cast failure: p_operator is not the correct DeltaSpin type.
- *   This happens if set_operator() was not called with the correct type.
- *   Solution: Ensure set_operator() is called before cal_mi_lcao().
+ * The actual LCAO magnetic-moment implementation lives in
+ * deltaspin_lcao_mi.cpp as free functions over ScState; the member
+ * functions below only adapt the singleton's stored pointers.
  */
 
-/**
- * @brief Calculate atomic magnetic moments using real-space projection (LCAO basis).
- *
- * @details The DeltaSpin operator computes magnetic moments by projecting the
- * density matrix onto atomic orbitals. For each constrained atom:
- *   M_i = Tr[P_at * (rho_up - rho_dn)]  (nspin=2)
- *   M_i = Tr[P_at * rho_spinor]          (nspin=4, decomposed via Pauli matrices)
- *
- * @param step Current SCF iteration number (for logging)
- * @param print Whether to print moments (unused in this implementation)
- */
+#include "spin_constrain.h"
+
+#include "deltaspin_lcao_mi.h"
+#include "source_lcao/module_operator_lcao/dspin_lcao.h"
+#include "source_estate/module_dm/density_matrix.h"
+
 template <>
 void spinconstrain::SpinConstrain<std::complex<double>>::cal_mi_lcao(const int& step, bool print)
 {
-    ModuleBase::TITLE("module_deltaspin", "cal_mi_lcao");
-    ModuleBase::timer::start("spinconstrain::SpinConstrain", "cal_mi_lcao");
-    // Reset Mi before calculation
-    this->zero_Mi();
-    const hamilt::HContainer<double>* dmr = this->dm_->get_DMR_pointer(1);
-    std::vector<double> moments;
-    if(this->nspin_==2)
-    {
-        // Switch to spin-difference density matrix (rho_up - rho_dn)
-        this->dm_->switch_dmr(2);
-
-        // Compute moments via DeltaSpin operator
-        moments = static_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, double>>*>(this->p_operator)->cal_moment(dmr, this->get_constrain());
-
-        // Switch back to total density matrix
-        this->dm_->switch_dmr(0);
-
-        // For nspin=2, only z-component is meaningful
-        for(int iat=0;iat<this->Mi_.size();iat++)
-        {
-            this->Mi_[iat].x = 0.0;
-            this->Mi_[iat].y = 0.0;
-            this->Mi_[iat].z = moments[iat];
-        }
-    }
-    else if(this->nspin_==4)
-    {
-        // For nspin=4, moments array contains interleaved [Mx, My, Mz] per atom
-        moments = static_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>*>(this->p_operator)->cal_moment(dmr, this->get_constrain());
-        for(int iat=0;iat<this->Mi_.size();iat++)
-        {
-            this->Mi_[iat].x = moments[iat*3];
-            this->Mi_[iat].y = moments[iat*3+1];
-            this->Mi_[iat].z = moments[iat*3+2];
-        }
-    }
-
-    ModuleBase::timer::end("spinconstrain::SpinConstrain", "cal_mi_lcao");
+    lcao::cal_mi_lcao(this->state_, this->p_operator, this->dm_, step, print);
 }
 
-#endif
+// cal_mi_lcao<double> stub lives in template_helpers.cpp (single definition).
 
-// cal_mi_pw() has been moved to source/source_pw/module_pwdft/deltaspin_pw_impl.cpp
-// because it depends on PW-specific OnsiteProjector.
-
-/// @brief Set the DeltaSpin operator pointer for LCAO magnetic moment calculation
 template <>
 void spinconstrain::SpinConstrain<std::complex<double>>::set_operator(
     hamilt::Operator<std::complex<double>>* op_in)
 {
     this->p_operator = op_in;
 }
 
-/// @brief Set the DeltaSpin operator pointer (double specialization)
 template <>
 void spinconstrain::SpinConstrain<double>::set_operator(
     hamilt::Operator<double>* op_in)
 {
     this->p_operator = op_in;
 }
+
+#endif // __LCAO

source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp

@@ -4,6 +4,7 @@
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_io/module_parameter/parameter.h"
 #include "spin_constrain.h"
+#include "deltaspin_pw_mi.h"
 #include "mi_tools.h"
 #include "source_pw/module_pwdft/onsite_proj.h"
 #include "source_base/parallel_reduce.h"
@@ -105,19 +106,19 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
                                                             PARAM.inp.nbands,
                                                             PARAM.inp.nelec,
                                                             PARAM.inp.device == "gpu");
-        if (this->nspin_ == 2)
+        if (this->state_.nspin_ == 2)
         {
             dynamic_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, double>>*>(this->p_operator)
                 ->update_lambda();
         }
-        else if (this->nspin_ == 4)
+        else if (this->state_.nspin_ == 4)
         {
             dynamic_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>*>(
                 this->p_operator)
                 ->update_lambda();
         }
         // Diagonalization without updating charge density (last param = true means skip charge update)
-        hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->dm_, *this->pelec->charge, this->nspin_, true);
+        hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->dm_, *this->pelec->charge, this->state_.nspin_, true);
         elecstate::calculate_weights(this->pelec->ekb,
                                      this->pelec->wg,
                                      this->pelec->klist,
@@ -164,23 +165,21 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
                 becp_tmp.resize(size_becp * nk);
                 std::vector<std::complex<double>> h_tmp(nbands * nbands), s_tmp(nbands * nbands);
                 int initial_hs = 0;
-                if(this->sub_h_save == nullptr)
+                if(!this->pw_cache_.allocated())
                 {
                     // FIRST CALL: save subspace data for reuse across lambda steps
                     initial_hs = 1;
-                    this->sub_h_save = new std::complex<double>[nbands * nbands * nk];
-                    this->sub_s_save = new std::complex<double>[nbands * nbands * nk];
-                    this->becp_save = new std::complex<double>[size_becp * nk];
-                    this->lambda_in_sub_ = this->lambda_;
+                    this->pw_cache_.allocate_cpu(nbands, nk, size_becp);
+                    this->pw_cache_.lambda_in_sub() = this->state_.lambda_;
                 }
                 for (int ik = 0; ik < nk; ++ik)
                 {
 
                     psi_t->fix_k(ik);
 
-                    std::complex<double>* h_k = this->sub_h_save + ik * nbands * nbands;
-                    std::complex<double>* s_k = this->sub_s_save + ik * nbands * nbands;
-                    std::complex<double>* becp_k = this->becp_save + ik * size_becp;
+                    std::complex<double>* h_k = this->pw_cache_.h_k(ik, nbands);
+                    std::complex<double>* s_k = this->pw_cache_.s_k(ik, nbands);
+                    std::complex<double>* becp_k = this->pw_cache_.becp_k(ik, size_becp);
                     if(initial_hs)
                     {
                         /// Compute H(k) and extract subspace matrices for this k-point
@@ -191,7 +190,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
                     memcpy(h_tmp.data(), h_k, sizeof(std::complex<double>) * nbands * nbands);
                     memcpy(s_tmp.data(), s_k, sizeof(std::complex<double>) * nbands * nbands);
                     // Apply DeltaSpin correction (skip for initialization step i_step=-1)
-                    if (i_step != -1) this->calculate_delta_hcc(h_tmp.data(), becp_k, this->lambda_.data(), nbands, nkb, nh_iat, ik, true);
+                    if (i_step != -1) pw::calculate_delta_hcc(this->state_, this->pw_cache_, this->pelec, h_tmp.data(), becp_k, this->state_.lambda_.data(), nbands, nkb, nh_iat, ik, true);
 
                     // Diagonalize in subspace, update becp (response wavefunctions)
                     hsolver::DiagoIterAssist<std::complex<double>>::diag_responce(h_tmp.data(),
@@ -224,23 +223,21 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
                 base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(h_tmp, nbands * nbands);
                 base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(s_tmp, nbands * nbands);
                 int initial_hs = 0;
-                if(this->sub_h_save == nullptr)
+                if(!this->pw_cache_.allocated())
                 {
                     initial_hs = 1;
-                    base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(this->sub_h_save, nbands * nbands * nk);
-                    base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(this->sub_s_save, nbands * nbands * nk);
-                    base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(this->becp_save, size_becp * nk);
-                    this->lambda_in_sub_ = this->lambda_;
+                    this->pw_cache_.allocate_gpu(nbands, nk, size_becp);
+                    this->pw_cache_.lambda_in_sub() = this->state_.lambda_;
                 }
                 std::complex<double>* becp_pointer = nullptr;
                 base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(becp_pointer, size_becp);
                 for (int ik = 0; ik < nk; ++ik)
                 {
                     psi_t->fix_k(ik);
 
-                    std::complex<double>* h_k = this->sub_h_save + ik * nbands * nbands;
-                    std::complex<double>* s_k = this->sub_s_save + ik * nbands * nbands;
-                    std::complex<double>* becp_k = this->becp_save + ik * size_becp;
+                    std::complex<double>* h_k = this->pw_cache_.h_k(ik, nbands);
+                    std::complex<double>* s_k = this->pw_cache_.s_k(ik, nbands);
+                    std::complex<double>* becp_k = this->pw_cache_.becp_k(ik, size_becp);
                     if(initial_hs)
                     {
                         hamilt_t->updateHk(ik);
@@ -249,7 +246,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
                     }
                     base_device::memory::synchronize_memory_op<std::complex<double>, base_device::DEVICE_GPU, base_device::DEVICE_GPU>()(h_tmp, h_k, nbands * nbands);
                     base_device::memory::synchronize_memory_op<std::complex<double>, base_device::DEVICE_GPU, base_device::DEVICE_GPU>()(s_tmp, s_k, nbands * nbands);
-                    if (i_step != -1) this->calculate_delta_hcc(h_tmp, becp_k, this->lambda_.data(), nbands, nkb, nh_iat, ik, true);
+                    if (i_step != -1) pw::calculate_delta_hcc(this->state_, this->pw_cache_, this->pelec, h_tmp, becp_k, this->state_.lambda_.data(), nbands, nkb, nh_iat, ik, true);
 
                     hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_GPU>::diag_responce(h_tmp,
                                                                                   s_tmp,
@@ -280,15 +277,15 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
             for (int ik = 0; ik < nk; ik++)
             {
                 const std::complex<double>* becp = &becp_tmp[ik * size_becp];
-                const int spin_sign = (this->npol_ == 2) ? 1 : this->get_spin_sign(ik);
-                accumulate_Mi_from_becp(becp, nkb, nbands, this->npol_, spin_sign,
-                    &this->pelec->wg(ik, 0), nh_iat, this->Mi_);
+                const int spin_sign = (this->state_.npol_ == 2) ? 1 : this->get_spin_sign(ik);
+                accumulate_Mi_from_becp(becp, nkb, nbands, this->state_.npol_, spin_sign,
+                    &this->pelec->wg(ik, 0), nh_iat, this->state_.Mi_);
             }
             // MPI reduction: sum Mi across all k-pool ranks
             Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar,
                                                     GlobalV::NPROC_IN_POOL,
-                                                    &(this->Mi_[0][0]),
-                                                    3 * this->Mi_.size());
+                                                    &(this->state_.Mi_[0][0]),
+                                                    3 * this->state_.Mi_.size());
         }
     }
     ModuleBase::timer::end("spinconstrain::SpinConstrain", "cal_mw_from_lambda");
@@ -328,12 +325,14 @@ void spinconstrain::SpinConstrain<std::complex<double>>::update_psi_charge(const
     {
         if (PARAM.inp.device == "cpu")
         {
-            this->update_psi_charge_pw_cpu(delta_lambda, pw_solve, full_update);
+            pw::update_psi_charge_pw_cpu(this->state_, this->pw_cache_, this->psi, this->p_hamilt,
+                                         this->pelec, this->pw_wfc_, delta_lambda, pw_solve, full_update);
         }
 #if ((defined __CUDA) || (defined __ROCM))
         else
         {
-            this->update_psi_charge_pw_gpu(delta_lambda, pw_solve, full_update);
+            pw::update_psi_charge_pw_gpu(this->state_, this->pw_cache_, this->psi, this->p_hamilt,
+                                         this->pelec, this->pw_wfc_, delta_lambda, pw_solve, full_update);
         }
 #endif
     }