Skip to content

refactor(deltaspin): encapsulate ScState lambda and Mi #8112

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

Merged
mohanchen merged 4 commits into from
Oct 10, 2026
Merged

refactor(deltaspin): encapsulate ScState lambda and Mi #8112

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
34456f5
refactor(deltaspin): encapsulate ScState::lambda_ and Mi_ as private …
hujieting Oct 10, 2026
519c755
refactor(deltaspin): rename mutable accessors per review
Oct 10, 2026
276087e
refactor(deltaspin): unify ScState accessor naming
Oct 10, 2026
6e479f7
merge: resolve develop conflicts for deltaspin encapsulation
hujieting Oct 10, 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
14 changes: 7 additions & 7 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 @@ -175,7 +175,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
// FIRST CALL: save subspace data for reuse across lambda steps
initial_hs = 1;
this->pw_cache_.allocate_cpu(nbands, nk, size_becp);
this->pw_cache_.lambda_in_sub() = this->state_.lambda_;
this->pw_cache_.lambda_in_sub() = this->state_.get_sc_lambda();
}
for (int ik = 0; ik < nk; ++ik)
{
Expand All @@ -199,7 +199,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) 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);
if (i_step != -1) pw::calculate_delta_hcc(this->state_, this->pw_cache_, this->pelec, h_tmp.data(), becp_k, this->state_.get_sc_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 @@ -237,7 +237,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
{
initial_hs = 1;
this->pw_cache_.allocate_gpu(nbands, nk, size_becp);
this->pw_cache_.lambda_in_sub() = this->state_.lambda_;
this->pw_cache_.lambda_in_sub() = this->state_.get_sc_lambda();
}
std::complex<double>* becp_pointer = nullptr;
base_device::memory::resize_memory_op<std::complex<double>, base_device::DEVICE_GPU>()(becp_pointer, size_becp);
Expand All @@ -261,7 +261,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) 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);
if (i_step != -1) pw::calculate_delta_hcc(this->state_, this->pw_cache_, this->pelec, h_tmp, becp_k, this->state_.get_sc_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 @@ -294,13 +294,13 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
const std::complex<double>* becp = &becp_tmp[ik * size_becp];
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_);
&this->pelec->wg(ik, 0), nh_iat, this->state_.get_mi());
}
// MPI reduction: sum Mi across all k-pool ranks
Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar,
GlobalV::NPROC_IN_POOL,
&(this->state_.Mi_[0][0]),
3 * this->state_.Mi_.size());
&(this->state_.get_mi()[0][0]),
3 * this->state_.get_mi().size());
}
}
ModuleBase::timer::end("spinconstrain::SpinConstrain", "cal_mw_from_lambda");
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/module_deltaspin/deltaspin_init.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ void init_sc_state(const ScInitParams& params, const UnitCell& ucell, ScState& s
// Step 4: Load target magnetic moments and initial lambda from UnitCell
// These are parsed from the STRU file's "sc_mag" and "lambda" keywords
state.set_target_mag(unitcell::get_target_mag(ucell.atoms, ucell.ntype, ucell.nat));
state.lambda_ = unitcell::get_lambda(ucell.atoms, ucell.ntype, ucell.nat);
state.get_lambda() = unitcell::get_lambda(ucell.atoms, ucell.ntype, ucell.nat);
state.constrain_ = unitcell::get_constrain(ucell.atoms, ucell.ntype, ucell.nat);

// Step 5: CRITICAL FIX for collinear spin (nspin=2)
Expand Down
28 changes: 14 additions & 14 deletions source/source_lcao/module_deltaspin/deltaspin_lcao_mi.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,22 +61,22 @@ void cal_mi_lcao(ScState& state,
dm->switch_dmr(0);

// For nspin=2, only z-component is meaningful
for (int iat = 0; iat < state.Mi_.size(); iat++)
for (int iat = 0; iat < state.get_mi().size(); iat++)
{
state.Mi_[iat].x = 0.0;
state.Mi_[iat].y = 0.0;
state.Mi_[iat].z = moments[iat];
state.get_mi()[iat].x = 0.0;
state.get_mi()[iat].y = 0.0;
state.get_mi()[iat].z = moments[iat];
}
}
else if (state.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>>>*>(p_operator)->cal_moment(dmr, state.get_constrain());
for (int iat = 0; iat < state.Mi_.size(); iat++)
for (int iat = 0; iat < state.get_mi().size(); iat++)
{
state.Mi_[iat].x = moments[iat * 3];
state.Mi_[iat].y = moments[iat * 3 + 1];
state.Mi_[iat].z = moments[iat * 3 + 2];
state.get_mi()[iat].x = moments[iat * 3];
state.get_mi()[iat].y = moments[iat * 3 + 1];
state.get_mi()[iat].z = moments[iat * 3 + 2];
}
}

Expand Down Expand Up @@ -167,15 +167,15 @@ void calculate_mw_from_orbitals(const std::vector<std::vector<std::vector<double
}
if (state.nspin_ == 2)
{
state.Mi_[iat].x = 0.0;
state.Mi_[iat].y = 0.0;
state.Mi_[iat].z = atom_mag;
state.get_mi()[iat].x = 0.0;
state.get_mi()[iat].y = 0.0;
state.get_mi()[iat].z = atom_mag;
}
else if (state.nspin_ == 4)
{
state.Mi_[iat].x = (std::abs(total_charge_soc[1]) < state.sc_thr_)? 0.0 : total_charge_soc[1];
state.Mi_[iat].y = (std::abs(total_charge_soc[2]) < state.sc_thr_)? 0.0 : total_charge_soc[2];
state.Mi_[iat].z = (std::abs(total_charge_soc[3]) < state.sc_thr_)? 0.0 : total_charge_soc[3];
state.get_mi()[iat].x = (std::abs(total_charge_soc[1]) < state.sc_thr_)? 0.0 : total_charge_soc[1];
state.get_mi()[iat].y = (std::abs(total_charge_soc[2]) < state.sc_thr_)? 0.0 : total_charge_soc[2];
state.get_mi()[iat].z = (std::abs(total_charge_soc[3]) < state.sc_thr_)? 0.0 : total_charge_soc[3];
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/module_deltaspin/deltaspin_lcao_mi.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ namespace lcao
* @details Uses the DeltaSpin operator to compute Tr(rho * mu) per atom.
* For nspin=2, extracts only the z-component. For nspin=4, extracts
* all three components from the interleaved 4-component spinor density matrix.
* Results are stored in state.Mi_ (indexed by global atom index iat).
* Results are stored in the state's magnetic moments (indexed by global atom index iat).
*
* @param state Constraint state (Mi_ written, indexing maps read)
* @param p_operator Base pointer to DeltaSpin<OperatorLCAO<...>>; nullptr aborts
Expand Down
14 changes: 7 additions & 7 deletions source/source_lcao/module_deltaspin/deltaspin_pw_mi.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ void cal_mi_pw(ScState& state,
int nkb = onsite_p->get_tot_nproj();
const int spin_sign = (npol == 2) ? 1 : spin_sign_at(state, pelec, ik);
accumulate_Mi_from_becp(becp, nkb, nbands, npol, spin_sign,
&pelec->wg(ik, 0), &onsite_p->get_nh(0), state.Mi_);
&pelec->wg(ik, 0), &onsite_p->get_nh(0), state.get_mi());
}
}
#if ((defined __CUDA) || (defined __ROCM))
Expand All @@ -105,12 +105,12 @@ void cal_mi_pw(ScState& state,
int nkb = onsite_p->get_size_becp() / nbands / npol;
const int spin_sign = (npol == 2) ? 1 : spin_sign_at(state, pelec, ik);
accumulate_Mi_from_becp(becp, nkb, nbands, npol, spin_sign,
&pelec->wg(ik, 0), &onsite_p->get_nh(0), state.Mi_);
&pelec->wg(ik, 0), &onsite_p->get_nh(0), state.get_mi());
}
}
#endif
// MPI reduction: sum Mi across all k-pool ranks
Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, GlobalV::NPROC_IN_POOL, &(state.Mi_[0][0]), 3 * state.Mi_.size());
Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, GlobalV::NPROC_IN_POOL, &(state.get_mi()[0][0]), 3 * state.get_mi().size());

ModuleBase::timer::end("spinconstrain::SpinConstrain", "cal_mi_pw");
}
Expand Down Expand Up @@ -175,7 +175,7 @@ void calculate_delta_hcc(ScState& state,
// The coefficient array uses {up-up, down-up, up-down, down-down}.
// Then: ps_up = coeff0 * becp_up + coeff2 * becp_dn
// ps_dn = coeff1 * becp_up + coeff3 * becp_dn
for (size_t iat = 0; iat < state.Mi_.size(); iat++)
for (size_t iat = 0; iat < state.get_mi().size(); iat++)
{
const int nproj = nh_iat[iat];
const std::complex<double> coefficients0(effective_lambda[iat][2], 0.0);
Expand Down Expand Up @@ -206,7 +206,7 @@ void calculate_delta_hcc(ScState& state,
// ps = lambda_z * spin_sign * becp
// spin_sign = +1 for spin-up k-points, -1 for spin-down
const int spin_sign = spin_sign_at(state, pelec, ik);
for (size_t iat = 0; iat < state.Mi_.size(); iat++)
for (size_t iat = 0; iat < state.get_mi().size(); iat++)
{
const int nproj = nh_iat[iat];
double coefficients0 = effective_lambda[iat][2] * spin_sign;
Expand Down Expand Up @@ -320,7 +320,7 @@ void update_psi_charge_pw_cpu(ScState& state,
const ModuleBase::Vector3<double>* lambda_for_hcc = delta_lambda;
if (full_update)
{
lambda_for_hcc = state.lambda_.data();
lambda_for_hcc = state.get_sc_lambda().data();
}

// =============================================================
Expand Down Expand Up @@ -439,7 +439,7 @@ void update_psi_charge_pw_gpu(ScState& state,
const ModuleBase::Vector3<double>* lambda_for_hcc = delta_lambda;
if (full_update)
{
lambda_for_hcc = state.lambda_.data();
lambda_for_hcc = state.get_sc_lambda().data();
}

// STAGE 1: Subspace diagonalization for each k-point (GPU)
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/module_deltaspin/deltaspin_pw_mi.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ class SubspaceCache;
* via OnsiteProjector::overlap_proj_psi, then decompose becp into magnetic moments.
* Finally Mi is summed across all MPI k-pool ranks.
*
* @param state Constraint state; state.Mi_ is filled in place.
* @param state Constraint state; its magnetic moments are filled in place.
* @param psi PW wavefunctions (psi::Psi<std::complex<double>>*, passed as void*
* to keep this header free of the Device template parameter).
* @param pelec Electronic state (provides wg weights and k-list spin signs).
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/module_deltaspin/deltaspin_state.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -530,7 +530,7 @@ double ScState::get_current_sc_thr() const
}

/// get computed magnetic moments Mi per atom
const std::vector<ModuleBase::Vector3<double>>& ScState::get_Mi() const
const std::vector<ModuleBase::Vector3<double>>& ScState::get_mi() const
{
return this->Mi_;
}
Expand Down
16 changes: 13 additions & 3 deletions source/source_lcao/module_deltaspin/deltaspin_state.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,11 @@ class ScState
void set_constrain(const ModuleBase::Vector3<int>* constrain_in, int nat_in);
/// get sc_lambda
const std::vector<ModuleBase::Vector3<double>>& get_sc_lambda() const;
/// get mutable sc_lambda
std::vector<ModuleBase::Vector3<double>>& get_lambda()
{
return lambda_;
}
/// get target_mag
const std::vector<ModuleBase::Vector3<double>>& get_target_mag() const;
/// get constrain
Expand Down Expand Up @@ -144,7 +149,12 @@ class ScState
/// get sc_drop_thr
double get_sc_drop_thr() const;
/// get computed magnetic moments Mi per atom
const std::vector<ModuleBase::Vector3<double>>& get_Mi() const;
const std::vector<ModuleBase::Vector3<double>>& get_mi() const;
/// get mutable computed magnetic moments Mi per atom
std::vector<ModuleBase::Vector3<double>>& get_mi()
{
return Mi_;
}
/// get human-readable atom labels ("Fe_0", "Fe_1", ...) for table printing
const std::vector<std::string>& get_atomLabels() const;
/// Total number of orbitals across all constrained atoms
Expand All @@ -171,9 +181,7 @@ class ScState
* later refactoring step.
* =============================================================
*/
std::vector<ModuleBase::Vector3<double>> lambda_; ///< Lagrange multipliers (Ry/uB) per atom, 3 components
std::vector<ModuleBase::Vector3<double>> target_mag_; ///< Target magnetic moments (uB) per atom
std::vector<ModuleBase::Vector3<double>> Mi_; ///< Current computed magnetic moments (uB) per atom
std::vector<ModuleBase::Vector3<int>> constrain_; ///< Per-atom/component constraint flags: 0=free, 1=constrained
std::vector<std::string> atomLabels_; ///< Human-readable labels: "Fe_0", "Fe_1", etc.
int nspin_ = 0; ///< Spin type: 2=collinear, 4=non-collinear
Expand All @@ -191,6 +199,8 @@ class ScState
double tpiba = 0.0; ///< 2*pi/a lattice constant scaling factor, saved from UnitCell

private:
std::vector<ModuleBase::Vector3<double>> lambda_; ///< Lagrange multipliers (Ry/uB) per atom, 3 components
std::vector<ModuleBase::Vector3<double>> Mi_; ///< Current computed magnetic moments (uB) per atom
std::map<int, std::vector<ScAtomData>> ScData; ///< Raw constraint data indexed by element type (itype)
std::map<int, double> ScDecayGrad; ///< Gradient decay thresholds (uB^2/eV) per element type
std::vector<double> decay_grad_; ///< Gradient decay thresholds converted to uB^2/Ry, per element type
Expand Down
28 changes: 14 additions & 14 deletions source/source_lcao/module_deltaspin/lambda_loop.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,10 +122,10 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
// Compute initial magnetic moments and save starting state
// =============================================================
this->cal_mw_from_lambda(i_step);
spin = this->state_.Mi_;
spin = this->state_.get_mi();

// Save initial lambda: for unconstrained components (constrain==0), set to 0
where_fill_scalar_else_2d(this->state_.constrain_, 0, zero, this->state_.lambda_, initial_lambda);
where_fill_scalar_else_2d(this->state_.constrain_, 0, zero, this->state_.get_sc_lambda(), initial_lambda);

print_2d(" initial lambda (eV/uB): ", initial_lambda, this->state_.nspin_, ModuleBase::Ry_to_eV, ofs_running);
print_2d(" initial spin (uB): ", spin, this->state_.nspin_, 1.0, ofs_running);
Expand All @@ -143,7 +143,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
where_fill_scalar_2d(this->state_.constrain_, 0, zero, delta_lambda);

// lambda = initial_lambda + delta_lambda
add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->state_.lambda_);
add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->state_.get_lambda());

// [direction_only mode] Project out parallel component of lambda
// This keeps |lambda| -> 0, only constraining spin direction
Expand All @@ -155,26 +155,26 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out

if (norm > 1e-8) {
const ModuleBase::Vector3<double> dir = target / norm;
double parallel = this->state_.lambda_[ia].x*dir.x +
this->state_.lambda_[ia].y*dir.y +
this->state_.lambda_[ia].z*dir.z;
this->state_.lambda_[ia].x -= parallel * dir.x;
this->state_.lambda_[ia].y -= parallel * dir.y;
this->state_.lambda_[ia].z -= parallel * dir.z;
double parallel = this->state_.get_lambda()[ia].x*dir.x +
this->state_.get_lambda()[ia].y*dir.y +
this->state_.get_lambda()[ia].z*dir.z;
this->state_.get_lambda()[ia].x -= parallel * dir.x;
this->state_.get_lambda()[ia].y -= parallel * dir.y;
this->state_.get_lambda()[ia].z -= parallel * dir.z;
}
}

// Apply lambda and compute new magnetic moments
this->cal_mw_from_lambda(i_step, delta_lambda.data());
new_spin = this->state_.Mi_;
new_spin = this->state_.get_mi();

// Check if gradient dM/dlambda has decayed below threshold
bool GradLessThanBound = check_gradient_decay(*this, new_spin, spin, delta_lambda, dnu_last_step, false, ofs_running);
if (i_step >= this->state_.nsc_min_ && GradLessThanBound)
{
// Gradient has decayed: further optimization yields diminishing returns
// Apply the last successful step and exit
add_scalar_multiply_2d(initial_lambda, dnu_last_step, one, this->state_.lambda_);
add_scalar_multiply_2d(initial_lambda, dnu_last_step, one, this->state_.get_lambda());
this->update_psi_charge(dnu_last_step.data(), true, true);
#ifdef __MPI
duration = (double)(MPI_Wtime() - iterstart);
Expand Down Expand Up @@ -269,7 +269,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
if(PARAM.inp.basis_type == "pw")
{
pw::cal_mi_pw(this->state_, this->psi, this->pelec);
subtract_2d(this->state_.Mi_, this->state_.target_mag_, delta_spin);
subtract_2d(this->state_.get_mi(), this->state_.target_mag_, delta_spin);
where_fill_scalar_2d(this->state_.constrain_, 0, zero, delta_spin);
search = delta_spin;
for (int ia = 0; ia < nat; ia++)
Expand Down Expand Up @@ -344,13 +344,13 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
// Mask unconstrained components
where_fill_scalar_else_2d(this->state_.constrain_, 0, zero, delta_lambda, delta_lambda);
// Update lambda
add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->state_.lambda_);
add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->state_.get_lambda());

// =============================================================
// TRIAL STEP: compute Mi at trial position
// =============================================================
this->cal_mw_from_lambda(i_step, delta_lambda.data());
spin_plus = this->state_.Mi_;
spin_plus = this->state_.get_mi();

// Find optimal step size via linear interpolation
alpha_opt = cal_alpha_opt(*this, spin, spin_plus, alpha_trial);
Expand Down
Loading
Loading