From e2a8728f04841c8b195cc8ff7c89910d25fe4859 Mon Sep 17 00:00:00 2001 From: hujieting <3114847804@qq.com> Date: Fri, 30 Oct 2026 20:12:50 +0800 Subject: [PATCH] fix(mulliken): correct noncollinear My sign Correct the nspin=4 Mulliken My sign for the conj-first density matrix convention. Add a regression test covering nonzero My for several noncollinear spin directions. Keep the nspin=1/2 paths unchanged. --- .../module_mulliken/output_mulliken.cpp | 110 ++++++++++-------- .../module_mulliken/output_mulliken.h | 17 +++ .../unittests/test_output_mulliken.cpp | 108 +++++++++++++++++ 3 files changed, 189 insertions(+), 46 deletions(-) diff --git a/source/source_io/module_mulliken/output_mulliken.cpp b/source/source_io/module_mulliken/output_mulliken.cpp index d477c0fa1e5..dc9a4d108ec 100644 --- a/source/source_io/module_mulliken/output_mulliken.cpp +++ b/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 Output_Mulliken::Output_Mulliken(Output_Sk* output_sk, Output_DMK* output_dmk, @@ -423,52 +486,7 @@ void Output_Mulliken::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_); } } diff --git a/source/source_io/module_mulliken/output_mulliken.h b/source/source_io/module_mulliken/output_mulliken.h index 933d721aa9e..3c8aaf3d65d 100644 --- a/source/source_io/module_mulliken/output_mulliken.h +++ b/source/source_io/module_mulliken/output_mulliken.h @@ -17,6 +17,23 @@ class Output_Sk; template 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 class Output_Mulliken diff --git a/source/source_io/unittests/test_output_mulliken.cpp b/source/source_io/unittests/test_output_mulliken.cpp index 2d4b0800d02..854211f42ea 100644 --- a/source/source_io/unittests/test_output_mulliken.cpp +++ b/source/source_io/unittests/test_output_mulliken.cpp @@ -8,6 +8,10 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include +#include +#include + template class OutputMullikenTest : public testing::Test { @@ -101,6 +105,110 @@ TYPED_TEST(OutputMullikenTest, nspin4) remove("./mulliken.txt"); } +namespace +{ +using cd = std::complex; + +/// 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); +} + +/// 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) {