Skip to content

Refactor DeltaSpin step 2: add Tests and refactor #7883

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mohanchen wants to merge 7 commits into
base: develop
Choose a base branch
from
Open

Refactor DeltaSpin step 2: add Tests and refactor #7883

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
3a96cd7
Test: add DeltaSpin integration tests with force/stress (PW + LCAO, n…
Aug 30, 2026
c4e9da6
Refactor DeltaSpin, steps 1-2: extract ScState and solver-independent…
Aug 30, 2026
6b52e0e
Refactor DeltaSpin, step 3: extract LCAO magnetic-moment path to free…
Aug 30, 2026
c49c67b
Fix CUDA build: qualify tpiba via state_ in update_psi_charge_pw_gpu
Aug 30, 2026
592e7d1
Refactor DeltaSpin, step 4: encapsulate PW subspace cache in Subspace…
Aug 30, 2026
dae8f04
Refactor DeltaSpin, step 5: move PW implementation back to module_del…
Aug 30, 2026
2178b2d
tests/03_NAO_multik/scf_deltaspin4: relax per-case threshold for know…
Aug 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions source/source_estate/test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,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
)
Expand Down
3 changes: 2 additions & 1 deletion source/source_io/module_ctrl/ctrl_output_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -214,7 +215,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);
}

Expand Down
6 changes: 4 additions & 2 deletions source/source_lcao/module_deltaspin/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
)

Expand Down
91 changes: 14 additions & 77 deletions source/source_lcao/module_deltaspin/cal_mw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
55 changes: 27 additions & 28 deletions source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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(),
Expand Down Expand Up @@ -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);
Expand All @@ -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,
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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
}
Expand Down
Loading
Loading