From 147228321c6061017c78080f234e6a13cab73fcf Mon Sep 17 00:00:00 2001 From: dyzheng Date: Wed, 16 Sep 2026 22:04:27 +0800 Subject: [PATCH 1/4] Fix(deltaspin): propagate current_spin through LCAO operator chain HamiltLCAO::updateHk() sets current_spin on the root operator via set_current_spin(isk[ik]), but the value was never forwarded to child operators. DeltaSpin did not toggle its own current_spin, so it always saw current_spin == 0 for nspin=2 and applied the same +lambda_z coefficient to both spin channels instead of +lambda_z/-lambda_z. The constraint therefore acted as a spin-independent potential and produced wrong magnetic moments and total energies. Propagate current_spin to the next operator in OperatorLCAO::init() before processing the current node, so every node in the chain shares the spin state set by the k-point loop. Regenerate tests/03_NAO_multik/scf_deltaspin2/result.ref, whose previous values encoded the buggy result. --- .../source_lcao/module_operator_lcao/operator_lcao.cpp | 8 ++++++++ tests/03_NAO_multik/scf_deltaspin2/result.ref | 9 +++++---- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp index f0560e69492..c1829b25d10 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/operator_lcao.cpp @@ -78,6 +78,14 @@ void OperatorLCAO::init(const int ik_in) { this->hR->set_zero(); } } + // propagate current_spin to next operator so all nodes in the chain + // share the same spin state, set by HamiltLCAO::updateHk via set_current_spin. + // This is done before processing this node so that children receive the + // correct spin regardless of any local toggling that may happen inside + // this operator's contributeHR(). + if (this->next_op != nullptr) { + dynamic_cast*>(this->next_op)->current_spin = this->current_spin; + } switch (this->cal_type) { case calculation_type::lcao_overlap: { // cal_type=lcao_overlap refer to overlap matrix operators, which are diff --git a/tests/03_NAO_multik/scf_deltaspin2/result.ref b/tests/03_NAO_multik/scf_deltaspin2/result.ref index fc858235e90..ea1d55e48c8 100644 --- a/tests/03_NAO_multik/scf_deltaspin2/result.ref +++ b/tests/03_NAO_multik/scf_deltaspin2/result.ref @@ -1,4 +1,5 @@ -etotref -6762.435776188675 -etotperatomref -3381.217888094338 -totalforceref 63.230574 -totalstressref 2916.957427 +etotref -6789.4013179079911424 +etotperatomref -3394.7006589540 +totalforceref 67.888296 +totalstressref 6225.500296 +totaltimeref 2.43 From 622cd2bb3ccc7031207d1f8856cb3b9e5c0c5528 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Wed, 16 Sep 2026 22:08:38 +0800 Subject: [PATCH 2/4] Fix(pw): correct nspin=2 DeltaSpin occupation output and per-atom labels cal_occupations() read the becp layout with the npol=2 stride (ib*2*nkb) and ignored the spin channel for nspin=2 (npol=1), so the projected atomic magnetization printed by print_orb_chg() was wrong for nspin=2. Index by the psi npol and store the spin-up/down occupancy in the up-up/down-down Pauli blocks so that Charge = occ[0]+occ[3] and Mag(z) = occ[0]-occ[3] print correctly; nspin=1 keeps a zero magnetization; nspin=4 keeps the interleaved spinor layout. Also build per-atom labels (Fe1, Fe2, ...) for the Total Magnetism / Magnetic force tables in print_orb_chg(), print_Mi() and print_Mag_Force(), instead of passing the per-type label vector (size ntype) to tables with nat rows. get_iat() is made const so the label helpers can read it from a const SpinConstrain reference. --- .../module_deltaspin/lambda_loop_helper.cpp | 24 ++++++++-- .../module_deltaspin/spin_constrain.h | 2 +- .../module_pwdft/onsite_proj_overlap.cpp | 47 ++++++++++++++----- .../module_pwdft/onsite_proj_print.cpp | 7 +-- 4 files changed, 61 insertions(+), 19 deletions(-) diff --git a/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp b/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp index 307dd312103..0aa0afa3de4 100644 --- a/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp +++ b/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp @@ -353,6 +353,14 @@ void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) const int nspin = sc.get_nspin(); const auto& Mi = sc.get_Mi(); const auto& atomLabel = sc.get_atomLabels(); + std::vector atom_labels_iat(nat); + for (const auto& it : sc.get_atomCounts()) + { + for (int ia = 0; ia < it.second; ++ia) + { + atom_labels_iat[sc.get_iat(it.first, ia)] = atomLabel[it.first] + std::to_string(ia + 1); + } + } std::vector mag_x(nat, 0.0); std::vector mag_y(nat, 0.0); std::vector mag_z(nat, 0.0); @@ -369,7 +377,7 @@ void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) { mag_z[iat] = Mi[iat].z; } - table << atomLabel << mag_z; + table << atom_labels_iat << mag_z; ofs_running << table.str() << std::endl; } else if (nspin == 4) @@ -387,7 +395,7 @@ void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) mag_y[iat] = Mi[iat].y; mag_z[iat] = Mi[iat].z; } - table << atomLabel << mag_x << mag_y << mag_z; + table << atom_labels_iat << mag_x << mag_y << mag_z; ofs_running << table.str() << std::endl; } } @@ -406,6 +414,14 @@ void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) const int nspin = sc.get_nspin(); const auto& lambda = sc.get_sc_lambda(); const auto& atomLabel = sc.get_atomLabels(); + std::vector atom_labels_iat(nat); + for (const auto& it : sc.get_atomCounts()) + { + for (int ia = 0; ia < it.second; ++ia) + { + atom_labels_iat[sc.get_iat(it.first, ia)] = atomLabel[it.first] + std::to_string(ia + 1); + } + } std::vector mag_force_x(nat, 0.0); std::vector mag_force_y(nat, 0.0); std::vector mag_force_z(nat, 0.0); @@ -422,7 +438,7 @@ void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) { mag_force_z[iat] = lambda[iat].z * ModuleBase::Ry_to_eV; } - table << atomLabel << mag_force_z; + table << atom_labels_iat << mag_force_z; ofs_running << table.str() << std::endl; } else if (nspin == 4) @@ -440,7 +456,7 @@ void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) mag_force_y[iat] = lambda[iat].y * ModuleBase::Ry_to_eV; mag_force_z[iat] = lambda[iat].z * ModuleBase::Ry_to_eV; } - table << atomLabel << mag_force_x << mag_force_y << mag_force_z; + table << atom_labels_iat << mag_force_x << mag_force_y << mag_force_z; ofs_running << table.str() << std::endl; } } diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h index 1af2c116273..343cef7bbe3 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.h +++ b/source/source_lcao/module_deltaspin/spin_constrain.h @@ -361,7 +361,7 @@ class SpinConstrain /// check atomCounts void check_atomCounts() const { state_.check_atomCounts(); } /// get iat - int get_iat(int itype, int atom_index) { return state_.get_iat(itype, atom_index); } + int get_iat(int itype, int atom_index) const { return state_.get_iat(itype, atom_index); } /// set nspin void set_nspin(int nspin) { state_.set_nspin(nspin); } /// get nspin diff --git a/source/source_pw/module_pwdft/onsite_proj_overlap.cpp b/source/source_pw/module_pwdft/onsite_proj_overlap.cpp index 989a6092148..b23b6bcbd40 100644 --- a/source/source_pw/module_pwdft/onsite_proj_overlap.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_overlap.cpp @@ -151,6 +151,7 @@ void projectors::OnsiteProjector::cal_occupations( // loop over k-points to calculate Mi of \sum_{k,i,l,m} const int nbands = psi_in->get_nbands(); + const int npol = psi_in->get_npol(); for(int ik = 0; ik < psi_in->get_nk(); ik++) { psi_in->fix_k(ik); @@ -159,16 +160,19 @@ void projectors::OnsiteProjector::cal_occupations( this->tabulate_atomic(ik); } // std::cout << __FILE__ << ":" << __LINE__ << " nbands = " << nbands << std::endl; - this->overlap_proj_psi( - nbands * psi_in->get_npol(), - psi_in->get_pointer()); + this->overlap_proj_psi(nbands * npol, psi_in->get_pointer()); const std::complex* becp_p = this->get_h_becp(); // becp(nbands*npol , nkb) // mag = wg * \sum_{nh}becp * becp int nkb = this->tot_nproj; - //nkb = 18; - //std::cout << "at " << __FILE__ << ": " << __LINE__ << " output nbands: " << nbands << std::endl; - //std::cout << "at " << __FILE__ << ": " << __LINE__ << " output nkb: " << nkb << std::endl; + // nspin=2 (npol=1): the spin-up and spin-down channels are separate + // k-points. Store spin-up occupancy in the up-up Pauli block (occ[0]) + // and spin-down occupancy in the down-down block (occ[3]) so that + // print_orb_chg() yields: + // Charge = occ[0] + occ[3], Mag(z) = occ[0] - occ[3] + // nspin=1 (npol=1): no spin polarization, split the occupancy evenly + // between occ[0] and occ[3] so that the printed magnetization is zero. + // nspin=4 (npol=2): both spinor components are interleaved per band. for(int ib = 0;ib::cal_occupations( for(int ih = 0; ih < nh; ih++) { const int occ_index = (begin_ih + ih) * 4; - const int index = ib*2*nkb + begin_ih + ih; - occs[occ_index] += weight * conj(becp_p[index]) * becp_p[index]; - occs[occ_index + 1] += weight * conj(becp_p[index]) * becp_p[index + nkb]; - occs[occ_index + 2] += weight * conj(becp_p[index + nkb]) * becp_p[index]; - occs[occ_index + 3] += weight * conj(becp_p[index + nkb]) * becp_p[index + nkb]; + if (npol == 1) + { + const int index = ib * nkb + begin_ih + ih; + const double occ = weight * (conj(becp_p[index]) * becp_p[index]).real(); + if (PARAM.inp.nspin == 2 && this->isk_ && this->isk_[ik] == 1) + { + occs[occ_index + 3] += occ; + } + else if (PARAM.inp.nspin == 1) + { + occs[occ_index] += 0.5 * occ; + occs[occ_index + 3] += 0.5 * occ; + } + else + { + occs[occ_index] += occ; + } + } + else + { + const int index = ib * 2 * nkb + begin_ih + ih; + occs[occ_index] += weight * conj(becp_p[index]) * becp_p[index]; + occs[occ_index + 1] += weight * conj(becp_p[index]) * becp_p[index + nkb]; + occs[occ_index + 2] += weight * conj(becp_p[index + nkb]) * becp_p[index]; + occs[occ_index + 3] += weight * conj(becp_p[index + nkb]) * becp_p[index + nkb]; + } } begin_ih += nh; } diff --git a/source/source_pw/module_pwdft/onsite_proj_print.cpp b/source/source_pw/module_pwdft/onsite_proj_print.cpp index adb6e1246c0..a2d5fec6679 100644 --- a/source/source_pw/module_pwdft/onsite_proj_print.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_print.cpp @@ -32,12 +32,13 @@ void print_orb_chg( /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); // parameters for mag output int occ_index = 0; + std::vector atom_labels_iat(ucell->nat); for(int iat=0; iatnat; iat++) { const int it = ucell->iat2it[iat]; - std::string atom_label = atom_labels[it]; int ia = ucell->iat2ia[iat]; - GlobalV::ofs_running << FmtCore::format("%-20s", atom_label+std::to_string(ia+1)) << std::endl; + atom_labels_iat[iat] = atom_labels[it] + std::to_string(ia+1); + GlobalV::ofs_running << FmtCore::format("%-20s", atom_labels_iat[iat]) << std::endl; std::vector sum(4, 0.0); int current_l = 1; std::vector charge_mag(4, 0.0); @@ -72,7 +73,7 @@ void print_orb_chg( GlobalV::ofs_running << std::endl; // Print magnetism table - print_mag_table(atom_labels, mag_x, mag_y, mag_z); + print_mag_table(atom_labels_iat, mag_x, mag_y, mag_z); } void print_mag_table( From f1e360ddae482c6f01fad59ba912ae5544213e7d Mon Sep 17 00:00:00 2001 From: dyzheng Date: Wed, 16 Sep 2026 22:16:39 +0800 Subject: [PATCH 3/4] Test: refresh stale DeltaSpin PW reference values The result.ref of these five PW DeltaSpin cases predates the DeltaSpin PW rework (they were last written in #7382) and no longer matched the code: 12/18 were off by 2.6/4.6 eV, while 19/21/41 differed only at the 1e-4-1e-7 eV level. Regenerate all five with the current code so the 17_DS_DFTU suite passes again. 12_PW_DS_S2_Z now agrees with the dedicated 01_PW/scf_deltaspin2 reference (-6369.19826815 eV), confirming the new value is the intended one. --- tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref | 6 +++--- tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref | 6 +++--- tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref | 6 +++--- tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref | 6 +++--- tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref | 6 +++--- 5 files changed, 15 insertions(+), 15 deletions(-) diff --git a/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref b/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref index 9dd83d8e706..dd48e508258 100644 --- a/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref +++ b/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref @@ -1,3 +1,3 @@ -etotref -6366.569118260046 -etotperatomref -3183.2845591300 -totaltimeref 1.97 +etotref -6369.198268154857 +etotperatomref -3184.5991340774 +totaltimeref 4.25 diff --git a/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref b/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref index 654cb15b3ee..72fb5c43b87 100644 --- a/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref +++ b/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref @@ -1,3 +1,3 @@ -etotref -6355.9855588350255857 -etotperatomref -3177.9927794175 -totaltimeref 2.89 +etotref -6360.5555597255606699 +etotperatomref -3180.2777798628 +totaltimeref 5.86 diff --git a/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref b/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref index e40be3f01bb..c3b83965fd3 100644 --- a/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref +++ b/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref @@ -1,3 +1,3 @@ -etotref -6360.5554588729937677 -etotperatomref -3180.277729436497 -totaltimeref 1.0 +etotref -6360.5555339529537378 +etotperatomref -3180.2777669765 +totaltimeref 9.56 diff --git a/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref b/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref index b06c7ee03a6..ae5a5a853dd 100644 --- a/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref +++ b/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref @@ -1,3 +1,3 @@ -etotref -6360.5554655414534864 -etotperatomref -3180.2777327707267 -totaltimeref 1.0 +etotref -6360.5554079455796455 +etotperatomref -3180.2777039728 +totaltimeref 6.65 diff --git a/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref b/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref index 8d127ea7298..eb51c776c6e 100644 --- a/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref +++ b/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref @@ -1,3 +1,3 @@ -etotref -6369.198254647004 -etotperatomref -3184.599127323502 -totaltimeref 1.0 +etotref -6369.198255055176 +etotperatomref -3184.5991275276 +totaltimeref 9.00 From b98aba704ff1d8354ce22171c68debf4d3e340b3 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Wed, 16 Sep 2026 22:36:35 +0800 Subject: [PATCH 4/4] Refactor(pw): pass nspin explicitly to cal_occupations The governance checker blocks PRs that increase GlobalV/GlobalC/PARAM usage. cal_occupations() read PARAM.inp.nspin twice; pass it as an explicit argument from ctrl_scf_pw() (which already holds the parsed Input_para) instead. Also keep the print_orb_chg() header line using the existing atom_label variable so the GlobalV::ofs_running line stays untouched. --- source/source_io/module_ctrl/ctrl_output_pw.cpp | 3 ++- source/source_pw/module_pwdft/onsite_proj.h | 4 +++- .../source_pw/module_pwdft/onsite_proj_overlap.cpp | 13 ++++++++----- source/source_pw/module_pwdft/onsite_proj_print.cpp | 5 +++-- 4 files changed, 16 insertions(+), 9 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index 58c88461641..ac525789c3e 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -227,7 +227,8 @@ void ModuleIO::ctrl_scf_pw(const int istep, { // float type has not been implemented auto* onsite_p = projectors::OnsiteProjector::get_instance(); onsite_p->cal_occupations(reinterpret_cast, Device>*>(stp.template get_psi_t()), - pelec->wg); + pelec->wg, + inp.nspin); } ModuleBase::timer::end("ModuleIO", "ctrl_scf_pw"); diff --git a/source/source_pw/module_pwdft/onsite_proj.h b/source/source_pw/module_pwdft/onsite_proj.h index a4181e4b4bb..0ede3a6e8c0 100644 --- a/source/source_pw/module_pwdft/onsite_proj.h +++ b/source/source_pw/module_pwdft/onsite_proj.h @@ -75,7 +75,9 @@ namespace projectors const ModuleBase::matrix& ekb); /// @brief calculate and print the occupations of all lm orbitals - void cal_occupations(const psi::Psi, Device>* psi, const ModuleBase::matrix& wg_in); + void cal_occupations(const psi::Psi, Device>* psi, + const ModuleBase::matrix& wg_in, + const int nspin_in); int get_size_becp() const { return size_becp; } std::complex* get_becp() const { return becp; } diff --git a/source/source_pw/module_pwdft/onsite_proj_overlap.cpp b/source/source_pw/module_pwdft/onsite_proj_overlap.cpp index b23b6bcbd40..aa2e280a0e9 100644 --- a/source/source_pw/module_pwdft/onsite_proj_overlap.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_overlap.cpp @@ -143,7 +143,8 @@ void projectors::OnsiteProjector::overlap_proj_psi( template void projectors::OnsiteProjector::cal_occupations( const psi::Psi, Device>* psi_in, - const ModuleBase::matrix& wg_in) + const ModuleBase::matrix& wg_in, + const int nspin_in) { ModuleBase::timer::start("OnsiteProj", "cal_occupation"); this->tabulate_atomic(0); @@ -187,11 +188,11 @@ void projectors::OnsiteProjector::cal_occupations( { const int index = ib * nkb + begin_ih + ih; const double occ = weight * (conj(becp_p[index]) * becp_p[index]).real(); - if (PARAM.inp.nspin == 2 && this->isk_ && this->isk_[ik] == 1) + if (nspin_in == 2 && this->isk_ && this->isk_[ik] == 1) { occs[occ_index + 3] += occ; } - else if (PARAM.inp.nspin == 1) + else if (nspin_in == 1) { occs[occ_index] += 0.5 * occ; occs[occ_index + 3] += 0.5 * occ; @@ -238,7 +239,8 @@ void projectors::OnsiteProjector::overlap_proj_ template void projectors::OnsiteProjector::cal_occupations( const psi::Psi, base_device::DEVICE_CPU>*, - const ModuleBase::matrix&); + const ModuleBase::matrix&, + const int); #if ((defined __CUDA) || (defined __ROCM)) template @@ -252,5 +254,6 @@ void projectors::OnsiteProjector::overlap_proj_ template void projectors::OnsiteProjector::cal_occupations( const psi::Psi, base_device::DEVICE_GPU>*, - const ModuleBase::matrix&); + const ModuleBase::matrix&, + const int); #endif diff --git a/source/source_pw/module_pwdft/onsite_proj_print.cpp b/source/source_pw/module_pwdft/onsite_proj_print.cpp index a2d5fec6679..6e97d133a81 100644 --- a/source/source_pw/module_pwdft/onsite_proj_print.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_print.cpp @@ -36,9 +36,10 @@ void print_orb_chg( for(int iat=0; iatnat; iat++) { const int it = ucell->iat2it[iat]; + std::string atom_label = atom_labels[it]; int ia = ucell->iat2ia[iat]; - atom_labels_iat[iat] = atom_labels[it] + std::to_string(ia+1); - GlobalV::ofs_running << FmtCore::format("%-20s", atom_labels_iat[iat]) << std::endl; + atom_labels_iat[iat] = atom_label + std::to_string(ia+1); + GlobalV::ofs_running << FmtCore::format("%-20s", atom_label+std::to_string(ia+1)) << std::endl; std::vector sum(4, 0.0); int current_l = 1; std::vector charge_mag(4, 0.0);