deepmodeling · mohanchen · Sep 18, 2026 · Sep 16, 2026 · Sep 16, 2026 · Sep 16, 2026
diff --git 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<double, Device>::get_instance();
         onsite_p->cal_occupations(reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(stp.template get_psi_t<T, Device>()),
-                                  pelec->wg);
+                                  pelec->wg,
+                                  inp.nspin);
     }
 
     ModuleBase::timer::end("ModuleIO", "ctrl_scf_pw");

diff --git 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<TK>& 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<std::string> 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<double> mag_x(nat, 0.0);
     std::vector<double> mag_y(nat, 0.0);
     std::vector<double> mag_z(nat, 0.0);
@@ -369,7 +377,7 @@ void print_Mi(const SpinConstrain<TK>& 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<TK>& 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<TK>& 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<std::string> 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<double> mag_force_x(nat, 0.0);
     std::vector<double> mag_force_y(nat, 0.0);
     std::vector<double> mag_force_z(nat, 0.0);
@@ -422,7 +438,7 @@ void print_Mag_Force(const SpinConstrain<TK>& 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<TK>& 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
@@ -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_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp
@@ -78,6 +78,14 @@ void OperatorLCAO<TK, TR>::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<OperatorLCAO<TK, TR>*>(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/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<std::complex<T>, Device>* psi, const ModuleBase::matrix& wg_in);
+        void cal_occupations(const psi::Psi<std::complex<T>, Device>* psi,
+                             const ModuleBase::matrix& wg_in,
+                             const int nspin_in);
 
         int get_size_becp() const { return size_becp; }
         std::complex<double>* 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
@@ -143,14 +143,16 @@ void projectors::OnsiteProjector<T, Device>::overlap_proj_psi(
 template<typename T, typename Device>
 void projectors::OnsiteProjector<T, Device>::cal_occupations(
         const psi::Psi<std::complex<T>, 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);
     std::vector<std::complex<double>> occs(this->tot_nproj * 4, 0.0);
 
     // loop over k-points to calculate Mi of \sum_{k,i,l,m}<Psi_{k,i}|alpha_{l,m}><alpha_{l,m}|Psi_{k,i}>
     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 +161,19 @@ void projectors::OnsiteProjector<T, Device>::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<double>* 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<nbands;ib++)
         {
             const double weight = wg_in(ik, ib);
@@ -179,11 +184,32 @@ void projectors::OnsiteProjector<T, Device>::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 (nspin_in == 2 && this->isk_ && this->isk_[ik] == 1)
+                        {
+                            occs[occ_index + 3] += occ;
+                        }
+                        else if (nspin_in == 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;
             }
@@ -213,7 +239,8 @@ void projectors::OnsiteProjector<double, base_device::DEVICE_CPU>::overlap_proj_
 template
 void projectors::OnsiteProjector<double, base_device::DEVICE_CPU>::cal_occupations(
     const psi::Psi<std::complex<double>, base_device::DEVICE_CPU>*,
-    const ModuleBase::matrix&);
+    const ModuleBase::matrix&,
+    const int);
 
 #if ((defined __CUDA) || (defined __ROCM))
 template
@@ -227,5 +254,6 @@ void projectors::OnsiteProjector<double, base_device::DEVICE_GPU>::overlap_proj_
 template
 void projectors::OnsiteProjector<double, base_device::DEVICE_GPU>::cal_occupations(
     const psi::Psi<std::complex<double>, 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
@@ -32,11 +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<std::string> atom_labels_iat(ucell->nat);
     for(int iat=0; iat<ucell->nat; iat++)
     {
         const int it = ucell->iat2it[iat];
         std::string atom_label = atom_labels[it];
         int ia = ucell->iat2ia[iat];
+        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<double> sum(4, 0.0);
         int current_l = 1;
@@ -72,7 +74,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(

diff --git 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
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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -1,3 +1,3 @@
-etotref -6369.198254647004
-etotperatomref -3184.599127323502
-totaltimeref 1.0
+etotref -6369.198255055176
+etotperatomref -3184.5991275276
+totaltimeref 9.00