Fix: correct My sign in noncollinear Mulliken magnetism

source/source_io/module_mulliken/output_mulliken.cpp

@@ -12,6 +12,69 @@
 namespace ModuleIO
 {
 
+void collect_MW_nspin4(ModuleBase::matrix& MecMulP,
+                       const ModuleBase::ComplexMatrix& mud,
+                       int nw,
+                       const Parallel_Orbitals* pv)
+{
+    for (size_t i = 0; i < nw; ++i)
+    {
+        const int index = i % 2;
+        if (!index)
+        {
+            const int j = i / 2;
+            const int k1 = 2 * j;
+            const int k2 = 2 * j + 1;
+            if (pv->in_this_processor(k1, k1))
+            {
+                const int ir = pv->global2local_row(k1);
+                const int ic = pv->global2local_col(k1);
+                MecMulP(0, j) += mud(ic, ir).real();
+                MecMulP(3, j) += mud(ic, ir).real();
+            }
+            // mud is DM.S with the conj-first density matrix built by
+            // dm_from_psi (DM_{ab} = sum_n w_n conj(c_{n,a}) c_{n,b}), i.e.
+            // DM = conj(P) for the physical one-body density matrix P. The
+            // up-down spin block therefore satisfies mud_{up,down} =
+            // conj(Q_{up,down}) for the physical Mulliken spin block Q, and
+            // the Pauli decomposition (see module_dm/dm_tools.cpp
+            // xyz_to_updown and module_deltaspin/mi_tools.h pauli_to_moment)
+            // reads
+            //   Mx = Re(mud_{up,down} + mud_{down,up})
+            //   My = Im(mud_{up,down} - mud_{down,up})
+            //   Mz = Re(mud_{up,up} - mud_{down,down}).
+            // Historically My was accumulated as
+            // -Im(mud_{up,down}) + Im(mud_{down,up}) = -My, i.e. with the
+            // wrong sign; that pre-existing flip is fixed here (see PR #7832,
+            // which documented it and deferred the fix to a separate PR).
+            if (pv->in_this_processor(k1, k2))
+            {
+                const int ir = pv->global2local_row(k1);
+                const int ic = pv->global2local_col(k2);
+                // note that mud is column major
+                MecMulP(1, j) += mud(ic, ir).real();
+                // M_y = Im(mud_{up,down} - mud_{down,up})
+                MecMulP(2, j) += mud(ic, ir).imag();
+            }
+            if (pv->in_this_processor(k2, k1))
+            {
+                const int ir = pv->global2local_row(k2);
+                const int ic = pv->global2local_col(k1);
+                MecMulP(1, j) += mud(ic, ir).real();
+                // M_y = Im(mud_{up,down} - mud_{down,up})
+                MecMulP(2, j) -= mud(ic, ir).imag();
+            }
+            if (pv->in_this_processor(k2, k2))
+            {
+                const int ir = pv->global2local_row(k2);
+                const int ic = pv->global2local_col(k2);
+                MecMulP(0, j) += mud(ic, ir).real();
+                MecMulP(3, j) -= mud(ic, ir).real();
+            }
+        }
+    }
+}
+
 template <typename TK>
 Output_Mulliken<TK>::Output_Mulliken(Output_Sk<TK>* output_sk,
                                      Output_DMK<TK>* output_dmk,
@@ -423,52 +486,7 @@ void Output_Mulliken<TK>::collect_MW(ModuleBase::matrix& MecMulP, const ModuleBa
     }
     else if (this->nspin_ == 4)
     {
-        for (size_t i = 0; i < nw; ++i)
-        {
-            const int index = i % 2;
-            if (!index)
-            {
-                const int j = i / 2;
-                const int k1 = 2 * j;
-                const int k2 = 2 * j + 1;
-                if (this->ParaV_->in_this_processor(k1, k1))
-                {
-                    const int ir = this->ParaV_->global2local_row(k1);
-                    const int ic = this->ParaV_->global2local_col(k1);
-                    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 dm_from_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);
-                    const int ic = this->ParaV_->global2local_col(k2);
-                    // note that mud is column major
-                    MecMulP(1, j) += mud(ic, ir).real();
-                    // M_y = i(M_{up,down} - M_{down,up}) = -(M_{up,down} - M_{down,up}).imag()
-                    MecMulP(2, j) -= mud(ic, ir).imag();
-                }
-                if (this->ParaV_->in_this_processor(k2, k1))
-                {
-                    const int ir = this->ParaV_->global2local_row(k2);
-                    const int ic = this->ParaV_->global2local_col(k1);
-                    MecMulP(1, j) += mud(ic, ir).real();
-                    // M_y = i(M_{up,down} - M_{down,up}) = -(M_{up,down} - M_{down,up}).imag()
-                    MecMulP(2, j) += mud(ic, ir).imag();
-                }
-                if (this->ParaV_->in_this_processor(k2, k2))
-                {
-                    const int ir = this->ParaV_->global2local_row(k2);
-                    const int ic = this->ParaV_->global2local_col(k2);
-                    MecMulP(0, j) += mud(ic, ir).real();
-                    MecMulP(3, j) -= mud(ic, ir).real();
-                }
-            }
-        }
+        collect_MW_nspin4(MecMulP, mud, nw, this->ParaV_);
     }
 }
 

source/source_io/module_mulliken/output_mulliken.h

@@ -17,6 +17,23 @@ class Output_Sk;
 template <typename TK>
 class Output_DMK;
 
+/**
+ * @brief Accumulate the nspin=4 Mulliken charge and Mx/My/Mz components.
+ *
+ * @details Decomposes one k-point mud = DM*S^T into the Pauli components using
+ * the conj-first DM convention (DM = conj(P), see dm_from_psi.h):
+ * Mx = Re(mud_ud + mud_du), My = Im(mud_ud - mud_du), Mz = Re(mud_uu - mud_dd).
+ *
+ * @param MecMulP accumulators [4 x nw/2]: 0 = charge, 1/2/3 = Mx/My/Mz
+ * @param mud     DM*S^T stored in ScaLAPACK column-major layout
+ * @param nw      total number of spinor orbitals (2 * number of spatial orbitals)
+ * @param pv      orbital distribution mapping
+ */
+void collect_MW_nspin4(ModuleBase::matrix& MecMulP,
+                       const ModuleBase::ComplexMatrix& mud,
+                       int nw,
+                       const Parallel_Orbitals* pv);
+
 /// @brief the output interface to write the Mulliken population charges
 template <typename TK>
 class Output_Mulliken

source/source_io/unittests/test_output_mulliken.cpp

@@ -8,6 +8,10 @@
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 
+#include <algorithm>
+#include <cmath>
+#include <complex>
+
 template <typename TK>
 class OutputMullikenTest : public testing::Test
 {
@@ -101,6 +105,110 @@ TYPED_TEST(OutputMullikenTest, nspin4)
     remove("./mulliken.txt");
 }
 
+namespace
+{
+using cd = std::complex<double>;
+
+/// spinor of the +1 eigenstate of mhat.sigma for a unit vector mhat
+void spinor_from_direction(const double mhat[3], cd c[2])
+{
+    double cos_theta = mhat[2];
+    cos_theta = std::max(-1.0, cos_theta);
+    cos_theta = std::min(1.0, cos_theta);
+    const double theta = std::acos(cos_theta);
+    const double phi = std::atan2(mhat[1], mhat[0]);
+    const double half_theta = 0.5 * theta;
+    const double cos_half = std::cos(half_theta);
+    const double sin_half = std::sin(half_theta);
+    const double cos_phi = std::cos(phi);
+    const double sin_phi = std::sin(phi);
+    c[0] = cd(cos_half, 0.0);
+    c[1] = sin_half * cd(cos_phi, sin_phi);
+}
+
+/// <c|sigma_i|c> with the standard Pauli matrices (sigma_y = [[0,-i],[i,0]])
+double pauli_expectation(const cd c[2], const int component)
+{
+    const cd sigma[3][2][2] = {
+        {{cd(0.0, 0.0), cd(1.0, 0.0)}, {cd(1.0, 0.0), cd(0.0, 0.0)}},
+        {{cd(0.0, 0.0), cd(0.0, -1.0)}, {cd(0.0, 1.0), cd(0.0, 0.0)}},
+        {{cd(1.0, 0.0), cd(0.0, 0.0)}, {cd(0.0, 0.0), cd(-1.0, 0.0)}},
+    };
+    cd value(0.0, 0.0);
+    for (int a = 0; a < 2; ++a)
+    {
+        for (int b = 0; b < 2; ++b)
+        {
+            value += std::conj(c[a]) * sigma[component][a][b] * c[b];
+        }
+    }
+    return value.real();
+}
+} // namespace
+
+// Regression test for the nspin=4 Mulliken M_y sign.
+// The stored DM is conj-first (DM_{ab} = conj(c_a) c_b, see dm_from_psi.h), so
+// the physical Pauli components of mud = DM.S (S = 1 here) are
+//   Mx = Re(mud_{up,down} + mud_{down,up})
+//   My = Im(mud_{up,down} - mud_{down,up})
+//   Mz = Re(mud_{up,up} - mud_{down,down}).
+// The pre-fix code accumulated -Im(mud_{up,down}) + Im(mud_{down,up}) = -My,
+// so this test fails before the fix and passes after it.
+TEST(OutputMullikenSpinTest, Nspin4MyFollowsPauliConvention)
+{
+    const double dirs[4][3] = {
+        {0.0, 1.0, 0.0},    // pure +y: transverse channel, maximal |My|
+        {0.0, -1.0, 0.0},   // pure -y
+        {0.6, 0.8, 0.0},    // in-plane moment with My > 0
+        {0.36, 0.48, -0.8}, // general three-dimensional direction
+    };
+
+    for (const auto& mhat : dirs)
+    {
+        cd c[2];
+        spinor_from_direction(mhat, c);
+
+        // DM in the conj-first storage convention of dm_from_psi
+        cd dm[4];
+        dm[0] = std::conj(c[0]) * c[0];
+        dm[1] = std::conj(c[0]) * c[1];
+        dm[2] = std::conj(c[1]) * c[0];
+        dm[3] = std::conj(c[1]) * c[1];
+
+        // one spatial orbital, one k-point, overlap S = 1: mud = DM * S^T = DM.
+        // mud is a ScaLAPACK column-major buffer; collect_MW_nspin4 accesses
+        // mud(ic, ir) as the (ir, ic) element of the logical matrix.
+        const int nw = 2;
+        Parallel_Orbitals pv;
+        pv.set_serial(nw, nw);
+        ModuleBase::ComplexMatrix mud(pv.ncol, pv.nrow, true);
+        for (int ir = 0; ir < nw; ++ir)
+        {
+            for (int ic = 0; ic < nw; ++ic)
+            {
+                const int idx = ir * 2 + ic;
+                mud(ic, ir) = dm[idx];
+            }
+        }
+
+        ModuleBase::matrix MecMulP(4, 1, true);
+        ModuleIO::collect_MW_nspin4(MecMulP, mud, nw, &pv);
+
+        const double expected_charge = 1.0;
+        const double expected_mx = pauli_expectation(c, 0);
+        const double expected_my = pauli_expectation(c, 1);
+        const double expected_mz = pauli_expectation(c, 2);
+        EXPECT_NEAR(MecMulP(0, 0), expected_charge, 1e-12);
+        EXPECT_NEAR(MecMulP(1, 0), expected_mx, 1e-12)
+            << "Mx wrong for direction (" << mhat[0] << "," << mhat[1] << "," << mhat[2] << ")";
+        EXPECT_NEAR(MecMulP(2, 0), expected_my, 1e-12)
+            << "My sign/value wrong (conj-first DM) for direction (" << mhat[0] << "," << mhat[1] << ","
+            << mhat[2] << ")";
+        EXPECT_NEAR(MecMulP(3, 0), expected_mz, 1e-12)
+            << "Mz wrong for direction (" << mhat[0] << "," << mhat[1] << "," << mhat[2] << ")";
+    }
+}
+
 #include "mpi.h"
 int main(int argc, char** argv)
 {