Fix: revert #7664 $m_y$ convention

source/source_cell/module_symmetry/test/CMakeLists.txt

@@ -21,4 +21,9 @@ AddTest(
   TARGET MODULE_CELL_SYMMETRY_rho_soc
   LIBS parameter base ${math_libs} device symmetry
   SOURCES symm_rho_soc_test.cpp
+    ${ABACUS_SOURCE_DIR}/source_estate/module_dm/density_matrix.cpp
+    ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/base_matrix.cpp
+    ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp
+    ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp
+    ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp
 )

source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp

@@ -6,6 +6,7 @@
 #include "../symmetry.h"
 #include "../symm_rot_spin.h"
 #include "source_cell/unitcell.h"
+#include "source_estate/module_dm/density_matrix.h" // real func_xyz_to_updown
 
 /************************************************
  *  unit test of Symmetry::rhog_symmetry_nspin4
@@ -187,30 +188,37 @@ TEST(RhogSymmetrySoc, GroupInvariance)
 }
 
 // ---------------------------------------------------------------------------
-// Coupling test (nonzero m_y): the spin-density rotation W used for the grid
-// symmetrization (spin_so3) MUST agree with the SU(2) rotation of the spinor
-// density block followed by the Pauli decomposition convention that the rest of
-// the code uses (func_xyz_to_updown, #7664):
-//     rho_0 = Re(uu+dd), rho_x = Re(ud+du),
-//     rho_y = -Im(ud) + Im(du),  rho_z = Re(uu-dd).
-// This is the check that the self-referential GroupInvariance test above cannot
-// make (it uses the same wspin as oracle). A y-channel handedness mismatch
-// between spin_so3 and this sigma_y=[[0,-i],[i,0]] convention shows up here.
+// Coupling test (nonzero m_y): the spin-density rotation W=spin_so3 used by psymmg_soc for the
+// grid symmetrization MUST agree with the SU(2) rotation of the physical spinor state followed by
+// the REAL func_xyz_to_updown extraction (which reads the conj-first stored DM, DM=conj(P), and
+// uses the bare +Im(ud)-Im(du)). This test now calls the actual func_xyz_to_updown rather than a
+// local re-implementation, so the grid-rotation and DM-extraction conventions cannot drift apart
+// silently (it fails on the #7664 m_y flip). The self-referential GroupInvariance test above
+// cannot catch this because it uses the same wspin as its own oracle.
 // ---------------------------------------------------------------------------
 namespace
 {
 using cd = std::complex<double>;
-// spinor block D = r0*I + m.sigma  (sigma_y = [[0,-i],[i,0]]); layout {uu,ud,du,dd}
+// PHYSICAL spinor block P = r0*I + m.sigma  (sigma_y = [[0,-i],[i,0]]); layout {uu,ud,du,dd}
 ModuleSymmetry::SpinRotation::Su2 block_from_pauli(double r0, double mx, double my, double mz)
 {
     return {cd(r0 + mz, 0.0), cd(mx, -my), cd(mx, my), cd(r0 - mz, 0.0)};
 }
-// func_xyz_to_updown extraction (NEW / #7664 convention); factor of 2 vs. m is harmless.
-void pauli_from_block(const ModuleSymmetry::SpinRotation::Su2& D, double& rx, double& ry, double& rz)
+// The runtime stores the DM conj-first (DM = conj(P), cal_dm_psi); this is what func_xyz_to_updown
+// actually consumes. Given a physical block P, the stored block is its element-wise conjugate.
+ModuleSymmetry::SpinRotation::Su2 stored_dm_from_phys(const ModuleSymmetry::SpinRotation::Su2& P)
 {
-    rx = (D[1] + D[2]).real();     // Re(ud+du)
-    ry = -D[1].imag() + D[2].imag(); // -Im(ud)+Im(du)
-    rz = (D[0] - D[3]).real();     // Re(uu-dd)
+    return {std::conj(P[0]), std::conj(P[1]), std::conj(P[2]), std::conj(P[3])};
+}
+// call the REAL func_xyz_to_updown on a 2x2 stored-DM block; return (m_x, m_y, m_z)
+ModuleBase::Vector3<double> real_extract(const ModuleSymmetry::SpinRotation::Su2& Dstored)
+{
+    const cd tmp[4] = {Dstored[0], Dstored[1], Dstored[2], Dstored[3]}; // {uu,ud,du,dd}
+    const int col_size = 2;
+    const int step_trace[4] = {0, 1, col_size, col_size + 1};
+    double out[4] = {0.0, 0.0, 0.0, 0.0}; // rho0/x/y/z written at icol=0
+    elecstate::DensityMatrix_Tools::func_xyz_to_updown<double>(tmp, 0, step_trace, out);
+    return ModuleBase::Vector3<double>(out[step_trace[1]], out[step_trace[2]], out[step_trace[3]]);
 }
 } // namespace
 
@@ -233,19 +241,30 @@ TEST(RhogSymmetrySoc, SpinConventionCoupling)
         EXPECT_NEAR(Wgrid.e31, Wpauli.e31, TOL) << "g=" << g; EXPECT_NEAR(Wgrid.e32, Wpauli.e32, TOL) << "g=" << g;
         EXPECT_NEAR(Wgrid.e33, Wpauli.e33, TOL) << "g=" << g;
 
-        // (2) rotate the spinor block, extract Pauli comps (new convention), compare to Wgrid*m
+        // (2) End-to-end with the REAL func_xyz_to_updown, exactly the runtime data flow:
+        //   physical block P(m) --conj--> stored DM (conj-first) --func_xyz_to_updown--> grid m.
+        //   Rotate the PHYSICAL block by the spinor SU(2) U (U P U^dagger, i.e. the physical state
+        //   rotation), conj to the stored block, extract again -> m'. psymmg_soc rotates the grid
+        //   components with Wgrid=spin_so3, so we must have  m' == Wgrid * m.  This catches any
+        //   mismatch (e.g. the #7664 m_y flip) between func_xyz_to_updown and spin_so3.
         for (const auto& m : mtest)
         {
-            const ModuleSymmetry::SpinRotation::Su2 D = block_from_pauli(2.0, m[0], m[1], m[2]);
-            const ModuleSymmetry::SpinRotation::Su2 Dp = ModuleSymmetry::SpinRotation::rotate_spin_block(D, U);
-            double rx, ry, rz;
-            pauli_from_block(Dp, rx, ry, rz);
-            // Wgrid acts on the physical m; the block carries 2*m, so compare against 2*(Wgrid*m).
-            const ModuleBase::Vector3<double> mv(m[0], m[1], m[2]);
-            const ModuleBase::Vector3<double> mrot = Wgrid * mv;
-            EXPECT_NEAR(rx, 2.0 * mrot.x, TOL) << "g=" << g;
-            EXPECT_NEAR(ry, 2.0 * mrot.y, TOL) << "g=" << g << " (y-channel handedness)";
-            EXPECT_NEAR(rz, 2.0 * mrot.z, TOL) << "g=" << g;
+            const ModuleSymmetry::SpinRotation::Su2 P  = block_from_pauli(2.0, m[0], m[1], m[2]);
+            const ModuleSymmetry::SpinRotation::Su2 Pp = ModuleSymmetry::SpinRotation::rotate_spin_block(P, U);
+
+            const ModuleBase::Vector3<double> mF  = real_extract(stored_dm_from_phys(P));
+            const ModuleBase::Vector3<double> mFp = real_extract(stored_dm_from_phys(Pp));
+
+            // (2a) extraction recovers the physical magnetization (block carries 2*m)
+            EXPECT_NEAR(mF.x, 2.0 * m[0], TOL) << "g=" << g;
+            EXPECT_NEAR(mF.y, 2.0 * m[1], TOL) << "g=" << g << " (m_y extraction)";
+            EXPECT_NEAR(mF.z, 2.0 * m[2], TOL) << "g=" << g;
+
+            // (2b) grid rotation spin_so3 agrees with the SU(2) block rotation + real extraction
+            const ModuleBase::Vector3<double> mrot = Wgrid * mF;
+            EXPECT_NEAR(mFp.x, mrot.x, TOL) << "g=" << g;
+            EXPECT_NEAR(mFp.y, mrot.y, TOL) << "g=" << g << " (y-channel handedness)";
+            EXPECT_NEAR(mFp.z, mrot.z, TOL) << "g=" << g;
         }
     }
 }

source/source_estate/module_charge/symm_rhog.cpp

@@ -90,26 +90,16 @@ void Symmetry_rho::psymmg_soc(std::complex<double>* rhog_x, std::complex<double>
 		// index [0,nrotk) unitary, [nrotk, nrotk+nrotk_anti) the spatial parts of the
 		// antiunitary elements Theta*g -- same layout as density_sym_ops().
 		const int na = symm.magnetic_nspin4 ? symm.nrotk_anti : 0;
-		// PR 7664 flipped the spinor->Pauli rho^y sign (func_xyz_to_updown), so the grid
-		// magnetization rho^y = chr.rho[2] consumed by rhog_symmetry_nspin4 now uses the
-		// standard sigma_y=[[0,-i],[i,0]] convention. Relative to that convention the
-		// pseudovector rotation applied to (rho^x,rho^y,rho^z) must be conjugated by
-		// S=diag(1,-1,1) (i.e. the y-channel handedness is flipped); without this the
-		// density symmetrization is out of sync only in the transverse y-channel.
-		// spin_so3() itself is left untouched (it is still the physical pseudovector 
-		// rotation used for magnetic-group detection on the STRU moments). 
-		auto yflip = [](ModuleBase::Matrix3 W) {
-			W.e12 = -W.e12; W.e21 = -W.e21; W.e23 = -W.e23; W.e32 = -W.e32; return W; };
 		std::vector<ModuleBase::Matrix3> wspin(symm.nrotk + na);
 		for (int i = 0; i < symm.nrotk; ++i)
 		{
 			const ModuleBase::Matrix3 gmatc = ilatvec * symm.gmatrix[i] * latvec;
-			wspin[i] = yflip(ModuleSymmetry::SpinRotation::spin_so3(gmatc));
+			wspin[i] = ModuleSymmetry::SpinRotation::spin_so3(gmatc);
 		}
 		for (int j = 0; j < na; ++j)
 		{
 			const ModuleBase::Matrix3 gmatc = ilatvec * symm.gmatrix_anti[j] * latvec;
-			wspin[symm.nrotk + j] = yflip(ModuleSymmetry::SpinRotation::spin_so3(gmatc));
+			wspin[symm.nrotk + j] = ModuleSymmetry::SpinRotation::spin_so3(gmatc);
 		}
 		return wspin;
 	};

source/source_estate/module_dm/density_matrix.cpp

@@ -655,7 +655,12 @@ void DensityMatrix_Tools::func_xyz_to_updown<double>(const std::complex<double>
 {
     target_DMR_mat[icol + step_trace[0]] = tmp[0].real() + tmp[3].real();  // rho_0 = (rho_upup + rho_downdown).real()
     target_DMR_mat[icol + step_trace[1]] = tmp[1].real() + tmp[2].real();  // rho_x = (rho_updown + rho_downup).real()
-    target_DMR_mat[icol + step_trace[2]] = -tmp[1].imag() + tmp[2].imag(); // rho_y = -Im(rho_updown) + Im(rho_downup)
+    // rho_y: the stored DM block is the complex conjugate of the physical 1-RDM P (cal_dm_psi builds
+    // DM_{ab}=sum conj(c_a) c_b = conj(P), so tmp[1]=DM_{ud}=conj(P_{ud})). Extracting m_y from the
+    // CONJUGATED block therefore carries the opposite sign of the bare-textbook formula; m_x/m_z read
+    // Re() and are conjugation-invariant. Using the bare formula (PR #7664) sign-flips m_y and quenches
+    // in-plane non-collinear moments (e.g. Mn3Sn 120-deg AFM); see issue #7831.
+    target_DMR_mat[icol + step_trace[2]] = tmp[1].imag() - tmp[2].imag();  // rho_y = Im(P_updown) - Im(P_downup)
     target_DMR_mat[icol + step_trace[3]] = tmp[0].real() - tmp[3].real();  // rho_z = (rho_upup - rho_downdown).real()
 }
 
@@ -664,7 +669,8 @@ void DensityMatrix_Tools::func_xyz_to_updown<std::complex<double>>(const std::co
 {
     target_DMR_mat[icol + step_trace[0]] = tmp[0] + tmp[3];                                         // rho_0 = (rho_upup + rho_downdown)
     target_DMR_mat[icol + step_trace[1]] = tmp[1] + tmp[2];                                         // rho_x = (rho_updown + rho_downup)
-    target_DMR_mat[icol + step_trace[2]] = ModuleBase::IMAG_UNIT * (tmp[1] - tmp[2]); // rho_y = i*(rho_updown - rho_downup)
+    // rho_y sign accounts for the conjugated stored DM block (conj(P)); see the <double> specialization above.
+    target_DMR_mat[icol + step_trace[2]] = -ModuleBase::IMAG_UNIT * (tmp[1] - tmp[2]); // rho_y = -i*(rho_updown - rho_downup)
     target_DMR_mat[icol + step_trace[3]] = tmp[0] - tmp[3];                                         // rho_z = (rho_upup - rho_downdown)
 }
 

source/source_estate/module_dm/test/CMakeLists.txt

@@ -48,3 +48,13 @@ AddTest(
   ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp
   ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp
 )
+
+AddTest(
+  TARGET MODULE_ESTATE_dm_soc_magnetization_roundtrip_test
+  LIBS parameter  base device
+  SOURCES test_soc_magnetization_roundtrip.cpp ../density_matrix.cpp ../density_matrix_io.cpp tmp_mocks.cpp
+  ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/base_matrix.cpp
+  ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp
+  ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp
+  ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp
+)

source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp

@@ -0,0 +1,146 @@
+#include "gtest/gtest.h"
+#include "source_estate/module_dm/density_matrix.h"
+
+#include <complex>
+#include <cmath>
+
+/************************************************************************
+ * Regression test for the nspin=4 (non-collinear/SOC) magnetization
+ * round-trip through the density-matrix pipeline.
+ *
+ * Physical invariant (must hold regardless of internal sign conventions):
+ *   the magnetization <sigma> of the occupied one-electron state that is
+ *   encoded in the density matrix must be recovered, with the CORRECT SIGN
+ *   in ALL THREE cartesian components, by func_xyz_to_updown().
+ *
+ * Why this test exists (regression for the #7664 nspin=4 m_y sign flip):
+ *   ABACUS builds the k-space DM as  DM_{ab} = sum_n w_n conj(c_{n,a}) c_{n,b}
+ *   (cal_dm_psi.cpp: the conj() is applied to the FIRST index a). Hence the
+ *   stored DM block is the complex conjugate of the physical 1-RDM P:
+ *       DM_{up,dn} = conj(c_up) c_dn = conj(P_{up,dn}).
+ *   Since m_x, m_z read Re() (conjugation-invariant) but m_y reads Im(),
+ *   ONLY m_y is sensitive to this conjugation. func_xyz_to_updown() must be
+ *   consistent with that stored convention. PR #7664 set the m_y extraction
+ *   to the "bare" textbook formula (valid for P, not for conj(P)), which
+ *   flips m_y for in-plane moments and quenches non-collinear order
+ *   (e.g. Mn3Sn 120-degree AFM). This test pins m_y down.
+ *
+ * The helper build_DM_block_as_cal_dm_psi() MUST mirror cal_dm_psi.cpp. If
+ * that convention is ever changed (e.g. the "upstream" fix that makes the DM
+ * hold the physical P), update the helper in the SAME commit so this test
+ * keeps asserting the physical invariant.
+ ************************************************************************/
+
+namespace
+{
+using cd = std::complex<double>;
+
+// spinor of the occupied state with <sigma> = mhat (the +1 eigenstate of mhat.sigma)
+void spinor_from_direction(const double mhat[3], cd c[2])
+{
+    // |+n> = (cos(th/2), sin(th/2) e^{i ph});  n=(sin th cos ph, sin th sin ph, cos th)
+    const double th = std::acos(std::max(-1.0, std::min(1.0, mhat[2])));
+    const double ph = std::atan2(mhat[1], mhat[0]);
+    c[0] = cd(std::cos(0.5 * th), 0.0);
+    c[1] = std::sin(0.5 * th) * cd(std::cos(ph), std::sin(ph));
+}
+
+// Build the 4 spinor-block DM elements EXACTLY as cal_dm_psi.cpp stores them:
+//   DM_{a,b} = sum_occ w * conj(c_a) * c_b     (conj on the first index)
+// layout tmp = {uu, ud, du, dd}
+void build_DM_block_as_cal_dm_psi(const cd c[2], double w, cd tmp[4])
+{
+    tmp[0] = w * std::conj(c[0]) * c[0]; // uu
+    tmp[1] = w * std::conj(c[0]) * c[1]; // ud
+    tmp[2] = w * std::conj(c[1]) * c[0]; // du
+    tmp[3] = w * std::conj(c[1]) * c[1]; // dd
+}
+
+// physical magnetization of a normalized spinor: m_i = <c| sigma_i |c>
+void physical_m(const cd c[2], double m[3])
+{
+    m[0] = 2.0 * std::real(std::conj(c[0]) * c[1]);
+    m[1] = 2.0 * std::imag(std::conj(c[0]) * c[1]);
+    m[2] = std::norm(c[0]) - std::norm(c[1]);
+}
+} // namespace
+
+TEST(SocMagnetizationRoundtrip, ExtractRecoversPhysicalMagnetization)
+{
+    // several magnetization directions, all with a nonzero transverse (y) part
+    const double dirs[5][3] = {
+        {0.0, 1.0, 0.0},                                  // pure +y (the critical case)
+        {0.0, -1.0, 0.0},                                 // pure -y (like Mn3Sn atom-1)
+        {0.6, 0.8, 0.0},                                  // in-plane 120-deg-like
+        {0.36, 0.48, -0.8},                               // general 3D
+        {-0.5, 0.5, 0.70710678},                          // general 3D
+    };
+
+    // step_trace for a single 2x2 spinor block written contiguously as a 2x2 (col_size=2)
+    const int col_size = 2;
+    const int step_trace[4] = {0, 1, col_size, col_size + 1};
+
+    for (const auto& mhat : dirs)
+    {
+        cd c[2];
+        spinor_from_direction(mhat, c);
+
+        double m_ref[3];
+        physical_m(c, m_ref); // the TRUE magnetization encoded in the state
+
+        cd tmp[4];
+        build_DM_block_as_cal_dm_psi(c, 1.0, tmp);
+
+        // 2x2 output buffer (row-major), func writes rho0/x/y/z into step_trace slots at icol=0
+        double out[4] = {0, 0, 0, 0};
+        elecstate::DensityMatrix_Tools::func_xyz_to_updown<double>(tmp, 0, step_trace, out);
+
+        const double mx = out[step_trace[1]];
+        const double my = out[step_trace[2]];
+        const double mz = out[step_trace[3]];
+
+        EXPECT_NEAR(mx, m_ref[0], 1e-10) << "m_x wrong for dir (" << mhat[0] << "," << mhat[1] << "," << mhat[2] << ")";
+        EXPECT_NEAR(my, m_ref[1], 1e-10) << "m_y SIGN/VALUE wrong (transverse channel, #7664 regression) for dir ("
+                                         << mhat[0] << "," << mhat[1] << "," << mhat[2] << ")";
+        EXPECT_NEAR(mz, m_ref[2], 1e-10) << "m_z wrong for dir (" << mhat[0] << "," << mhat[1] << "," << mhat[2] << ")";
+    }
+}
+
+// Same invariant for the <complex> (multi-k) specialization, which is changed identically.
+// For a single occupied state the 2x2 block is Hermitian, so the extracted Pauli components come
+// out real and must equal the physical magnetization; the imaginary parts must vanish.
+TEST(SocMagnetizationRoundtrip, ComplexSpecializationRecoversPhysicalMagnetization)
+{
+    const double dirs[4][3] = {
+        {0.0, 1.0, 0.0}, {0.0, -1.0, 0.0}, {0.6, 0.8, 0.0}, {0.36, 0.48, -0.8},
+    };
+    const int col_size = 2;
+    const int step_trace[4] = {0, 1, col_size, col_size + 1};
+
+    for (const auto& mhat : dirs)
+    {
+        cd c[2];
+        spinor_from_direction(mhat, c);
+        double m_ref[3];
+        physical_m(c, m_ref);
+
+        cd tmp[4];
+        build_DM_block_as_cal_dm_psi(c, 1.0, tmp);
+
+        cd out[4] = {cd(0, 0), cd(0, 0), cd(0, 0), cd(0, 0)};
+        elecstate::DensityMatrix_Tools::func_xyz_to_updown<std::complex<double>>(tmp, 0, step_trace, out);
+
+        EXPECT_NEAR(out[step_trace[1]].real(), m_ref[0], 1e-10) << "m_x";
+        EXPECT_NEAR(out[step_trace[2]].real(), m_ref[1], 1e-10) << "m_y (complex specialization)";
+        EXPECT_NEAR(out[step_trace[3]].real(), m_ref[2], 1e-10) << "m_z";
+        EXPECT_NEAR(out[step_trace[1]].imag(), 0.0, 1e-10);
+        EXPECT_NEAR(out[step_trace[2]].imag(), 0.0, 1e-10);
+        EXPECT_NEAR(out[step_trace[3]].imag(), 0.0, 1e-10);
+    }
+}
+
+int main(int argc, char** argv)
+{
+    testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}

source/source_io/module_mulliken/output_mulliken.cpp

@@ -438,6 +438,11 @@ void Output_Mulliken<TK>::collect_MW(ModuleBase::matrix& MecMulP, const ModuleBa
                     MecMulP(0, j) += mud(ic, ir).real();
                     MecMulP(3, j) += mud(ic, ir).real();
                 }
+                // WARNING (pre-existing, predates #7664, since 3.7.0 commit a339356): this M_y sign
+                // is suspect. mud is DM.S with the same conj-first DM convention as cal_dm_psi, so the
+                // bare Im formula here likely yields -m_y. It must NOT be used as an oracle for the DM
+                // convention (it is probably why #7664's DM m_y flip went unnoticed). Fix + unit test
+                // should be a separate PR after verifying against a case with nonzero in-plane moment.
                 if (this->ParaV_->in_this_processor(k1, k2))
                 {
                     const int ir = this->ParaV_->global2local_row(k1);

source/source_lcao/module_deltaspin/spin_constrain.cpp

@@ -142,7 +142,7 @@ int SpinConstrain<TK>::get_spin_sign(int ik) const
  * where P_at = sum_{l,m} |alpha_{l,m}><alpha_{l,m}| is the atomic projector.
  *
  * The magnetic moment is extracted via Pauli matrix traces:
- *   Mx = Re(occ[1] + occ[2]), My = -Im(occ[1] - occ[2]), Mz = Re(occ[0] - occ[3])
+ *   Mx = Re(occ[1] + occ[2]), My = Im(occ[1] - occ[2]), Mz = Re(occ[0] - occ[3])  (occ conj-first)
  *
  * @par Algorithm (npol=1, collinear):
  * Only the z-component (spin projection) is computed: