From c80303afa47489e345804bd1eb3255feed48dd32 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 18 Aug 2026 02:14:58 -0400 Subject: [PATCH 1/6] revert #7664 m_y convention Fix comment Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- source/source_estate/module_charge/symm_rhog.cpp | 14 ++------------ source/source_estate/module_dm/density_matrix.cpp | 10 ++++++++-- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/source/source_estate/module_charge/symm_rhog.cpp b/source/source_estate/module_charge/symm_rhog.cpp index dd489c221b0..4327ee68603 100644 --- a/source/source_estate/module_charge/symm_rhog.cpp +++ b/source/source_estate/module_charge/symm_rhog.cpp @@ -90,26 +90,16 @@ void Symmetry_rho::psymmg_soc(std::complex* rhog_x, std::complex // 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 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; }; diff --git a/source/source_estate/module_dm/density_matrix.cpp b/source/source_estate/module_dm/density_matrix.cpp index 586010844bb..6f780e51bdd 100644 --- a/source/source_estate/module_dm/density_matrix.cpp +++ b/source/source_estate/module_dm/density_matrix.cpp @@ -655,7 +655,12 @@ void DensityMatrix_Tools::func_xyz_to_updown(const std::complex { 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>(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 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) } From abc360070c6f9bf07ac44e3dc3a4b39d1253fb68 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 18 Aug 2026 03:30:20 -0400 Subject: [PATCH 2/6] test: revert scf_*_spin4 result.ref to pre-#7664 values These nspin=4 LCAO references were updated by #7664 (1457e26fe) to match its m_y sign convention in func_xyz_to_updown, which sign-flips the transverse magnetization (the stored DM block is conj(P); m_y reads Im so it flips). Commit fa8492f97 reverts that convention and restores the correct pre-#7664 physics (verified: Mn3Sn recovers the 120-deg AFM ground state, lower energy). scf_u_spin4 (mag[111]+U): force 14.360->11.332, stress 4334.0->4697.8, etot -6789.142->-6789.282 (the #7664 bug raised the energy by ~0.14 eV). scf_out_dos_spin4 / scf_angle_spin4: tiny (1e-4) reverts, near-collinear. 099_PW_DJ_SO left unchanged (PW path, not affected by the LCAO fix). Co-Authored-By: Claude Opus 4.8 --- tests/03_NAO_multik/scf_angle_spin4/result.ref | 8 ++++---- tests/03_NAO_multik/scf_out_dos_spin4/result.ref | 6 +++--- tests/03_NAO_multik/scf_u_spin4/result.ref | 8 ++++---- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/03_NAO_multik/scf_angle_spin4/result.ref b/tests/03_NAO_multik/scf_angle_spin4/result.ref index 5bcf1d29508..e1656f8f8ce 100644 --- a/tests/03_NAO_multik/scf_angle_spin4/result.ref +++ b/tests/03_NAO_multik/scf_angle_spin4/result.ref @@ -1,5 +1,5 @@ -etotref -6267.4651944939805617 -etotperatomref -3133.7325972470 -totalforceref 0.000008 -totalstressref 3912.920542 +etotref -6267.4651896196382950 +etotperatomref -3133.7325948098 +totalforceref 0.000000 +totalstressref 3912.920437 totaltimeref 15.08 diff --git a/tests/03_NAO_multik/scf_out_dos_spin4/result.ref b/tests/03_NAO_multik/scf_out_dos_spin4/result.ref index 27634ae24d1..e54b51d3064 100644 --- a/tests/03_NAO_multik/scf_out_dos_spin4/result.ref +++ b/tests/03_NAO_multik/scf_out_dos_spin4/result.ref @@ -1,6 +1,6 @@ -etotref -1964.0663947982982336 +etotref -1964.0663947982770878 etotperatomref -982.0331973991 -totalforceref 0.162298 -totalstressref 1877.059021 +totalforceref 0.162158 +totalstressref 1877.059089 totaldosref 38 totaltimeref 16.23 diff --git a/tests/03_NAO_multik/scf_u_spin4/result.ref b/tests/03_NAO_multik/scf_u_spin4/result.ref index ac66570b979..bdf978fbbe7 100644 --- a/tests/03_NAO_multik/scf_u_spin4/result.ref +++ b/tests/03_NAO_multik/scf_u_spin4/result.ref @@ -1,5 +1,5 @@ -etotref -6789.1423886377124290 -etotperatomref -3394.5711943189 -totalforceref 14.359774 -totalstressref 4333.997 +etotref -6789.2816406266510967 +etotperatomref -3394.6408203133 +totalforceref 11.331534 +totalstressref 4697.832232 totaltimeref 9.71 From 08246fb53ed3da195963a30978eb88e3db024ae2 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 18 Aug 2026 03:33:00 -0400 Subject: [PATCH 3/6] test: add nspin=4 magnetization round-trip regression test Guards the func_xyz_to_updown m_y convention against the #7664 regression. Asserts the physical invariant: the magnetization encoded in a DM spinor block (built with cal_dm_psi's convention DM=conj(P)) is recovered with the correct sign in ALL THREE components by func_xyz_to_updown. The m_y channel (Im) is the sensitive one; m_x/m_z (Re) are conj-invariant. Fails on the #7664 convention (m_y sign-flipped), passes after the fix. Co-Authored-By: Claude Opus 4.8 --- .../module_dm/test/CMakeLists.txt | 10 ++ .../test/test_soc_magnetization_roundtrip.cpp | 113 ++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp diff --git a/source/source_estate/module_dm/test/CMakeLists.txt b/source/source_estate/module_dm/test/CMakeLists.txt index 02812bfa7db..14db1ee4f94 100644 --- a/source/source_estate/module_dm/test/CMakeLists.txt +++ b/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 +) diff --git a/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp b/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp new file mode 100644 index 00000000000..24e4f530343 --- /dev/null +++ b/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp @@ -0,0 +1,113 @@ +#include "gtest/gtest.h" +#include "source_estate/module_dm/density_matrix.h" + +#include +#include + +/************************************************************************ + * 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 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 (see ref/MSG/2026-08-#7664-...): + * 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; + +// spinor of the occupied state with = 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 = +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(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] << ")"; + } +} + +int main(int argc, char** argv) +{ + testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} From 7322465dd4be02d7bdd15203d01c308a854ef3ba Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 18 Aug 2026 09:31:57 -0400 Subject: [PATCH 4/6] fix(soc): reconcile DeltaSpin/tests with reverted nspin=4 m_y convention MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Follow-up to #7832 (revert #7664 m_y convention) addressing review feedback. DeltaSpin (was left in the #7664 y-convention while the DM was reverted): - dspin_lcao.cpp cal_coeff_lambda: revert to lambda_ud = lambda_x + i*lambda_y. pauli_to_moment already reads the bare +Im (via #7748) and measures the physical +m_y from the conj-first occupation; #7664 flipped both the measurement and this drive (y-mirror), #7748 reverted only the measurement, so the constraint loop was driving the y-mirror of the target. Restoring the pre-#7664 drive re-pairs it with the (correct) measurement. - spin_constrain.h/.cpp: fix stale comments (occ is conj-first; My=+Im). - deltaspin_core_test.cpp: local re-impl and GeneralCase pinned My=-0.4 (#7664); corrected to the bare +0.4 so it stops enshrining the wrong sign. Tests: - test_soc_magnetization_roundtrip: add the (multi-k) specialization case; drop the dangling ref/MSG doc path from the comment. - symm_rho_soc_test: relabel SpinConventionCoupling — it validates spin_so3 in the PHYSICAL-rho frame, not the runtime conj-first DM path; note it does not pin the DM convention. Mulliken: - output_mulliken.cpp: comment flagging that M_y (pre-existing, since 3.7.0) likely reports -m_y from the conj-first DM and must not be used as an oracle for the DM convention; real fix + unit test deferred to a separate PR. Not covered here (need integration runs / separate PR): regenerate & report the nspin=4 reference values from post-PR code; nspin=4 DeltaSpin regression with nonzero m_y target; in-plane non-collinear symmetry integration test; nspin=4 DMR restart compatibility with #7664-era files. Co-Authored-By: Claude Opus 4.8 --- .../test/symm_rho_soc_test.cpp | 21 +++-- .../test/test_soc_magnetization_roundtrip.cpp | 35 +++++++- .../module_mulliken/output_mulliken.cpp | 5 ++ .../module_deltaspin/spin_constrain.cpp | 86 +++++++++++++++++++ .../module_deltaspin/spin_constrain.h | 30 +++++++ .../test/deltaspin_core_test.cpp | 11 ++- .../module_operator_lcao/dspin_lcao.cpp | 11 ++- 7 files changed, 183 insertions(+), 16 deletions(-) diff --git a/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp b/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp index 4f9c686e5e9..9fcc9e458d9 100644 --- a/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp +++ b/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp @@ -188,24 +188,29 @@ 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): +// symmetrization (spin_so3) MUST agree with the SU(2) rotation of a spinor density +// block followed by the standard sigma_y=[[0,-i],[i,0]] Pauli decomposition: // 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. +// NOTE: this uses the PHYSICAL spin-density block (D_ud = mx - i*my), so the textbook +// -Im(ud)+Im(du) is correct here. This is NOT the frame the runtime uses: the actual +// stored DM is conj-first (DM=conj(P), cal_dm_psi), so func_xyz_to_updown / psymmg_soc +// consume conj(P) and use the BARE +Im(ud)-Im(du) (see #7832 / the DM round-trip test in +// source_estate/module_dm/test). This test is a spin_so3 sanity check in the physical +// frame; it does NOT validate the DM-path sign and must not be read as pinning the #7664 +// convention. TODO: exercise the real func_xyz_to_updown/psymmg_soc instead of this local +// re-implementation so the two conventions cannot drift apart silently. // --------------------------------------------------------------------------- namespace { using cd = std::complex; -// spinor block D = r0*I + m.sigma (sigma_y = [[0,-i],[i,0]]); layout {uu,ud,du,dd} +// PHYSICAL spinor block D = 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. +// textbook Pauli extraction from the PHYSICAL block (distinct from func_xyz_to_updown, which +// reads the conj-first stored DM); factor of 2 vs. m is harmless. void pauli_from_block(const ModuleSymmetry::SpinRotation::Su2& D, double& rx, double& ry, double& rz) { rx = (D[1] + D[2]).real(); // Re(ud+du) diff --git a/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp b/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp index 24e4f530343..419c0211746 100644 --- a/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp +++ b/source/source_estate/module_dm/test/test_soc_magnetization_roundtrip.cpp @@ -13,7 +13,7 @@ * 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 (see ref/MSG/2026-08-#7664-...): + * 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: @@ -106,6 +106,39 @@ TEST(SocMagnetizationRoundtrip, ExtractRecoversPhysicalMagnetization) } } +// Same invariant for the (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>(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); diff --git a/source/source_io/module_mulliken/output_mulliken.cpp b/source/source_io/module_mulliken/output_mulliken.cpp index 3594bf1d04f..a687ddaf241 100644 --- a/source/source_io/module_mulliken/output_mulliken.cpp +++ b/source/source_io/module_mulliken/output_mulliken.cpp @@ -438,6 +438,11 @@ void Output_Mulliken::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); diff --git a/source/source_lcao/module_deltaspin/spin_constrain.cpp b/source/source_lcao/module_deltaspin/spin_constrain.cpp index fd9d1d83f08..003b2e422fd 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.cpp +++ b/source/source_lcao/module_deltaspin/spin_constrain.cpp @@ -130,6 +130,92 @@ int SpinConstrain::get_spin_sign(int ik) const return (this->pelec->klist->isk[ik] == 0) ? 1 : -1; } +<<<<<<< HEAD +======= +/** + * @brief Accumulate magnetic moments from projector coefficients (becp) for one k-point. + * + * @par Algorithm (npol=2, non-collinear): + * For each atom, compute the 2x2 occupation matrix from becp coefficients: + * occ[0] = sum_ih becp_up^*(ih) * becp_up(ih) = + * occ[1] = sum_ih becp_up^*(ih) * becp_dn(ih) = + * occ[2] = sum_ih becp_dn^*(ih) * becp_up(ih) = + * occ[3] = sum_ih becp_dn^*(ih) * becp_dn(ih) = + * where P_at = sum_{l,m} |alpha_{l,m}> + * Mz += weight * occ * spin_sign + * where spin_sign = +1 for spin-up, -1 for spin-down. + * + * @param becp Projector coefficients, layout: [ib * npol * nkb + spin * nkb + ih] + * @param nkb Total number of projectors across all atoms + * @param nbands Number of bands (occupied + unoccupied in the subspace) + * @param npol Number of spinor components (1 for collinear, 2 for non-collinear) + * @param ik K-point index (used for spin_sign lookup in collinear mode) + * @param wg_ik Band occupation weights for this k-point (from Fermi-Dirac) + * @param nh_iat Array of projector counts per atom: nh_iat[iat] = nproj for atom iat + */ +template +void SpinConstrain::accumulate_Mi_from_becp(const std::complex* becp, + int nkb, + int nbands, + int npol, + int ik, + const double* wg_ik, + const int* nh_iat) +{ + if (npol == 2) + { + for (int ib = 0; ib < nbands; ib++) + { + const double weight = wg_ik[ib]; + int begin_ih = 0; + for (int iat = 0; iat < static_cast(this->Mi_.size()); iat++) + { + std::complex occ[4] = {ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO}; + const int nh = nh_iat[iat]; + for (int ih = 0; ih < nh; ih++) + { + const int index = ib * 2 * nkb + begin_ih + ih; + occ[0] += conj(becp[index]) * becp[index]; + occ[1] += conj(becp[index]) * becp[index + nkb]; + occ[2] += conj(becp[index + nkb]) * becp[index]; + occ[3] += conj(becp[index + nkb]) * becp[index + nkb]; + } + this->Mi_[iat] += pauli_to_moment(occ, weight); + begin_ih += nh; + } + } + } + else // npol == 1 + { + const int sign = this->get_spin_sign(ik); + for (int ib = 0; ib < nbands; ib++) + { + const double weight = wg_ik[ib]; + int begin_ih = 0; + for (int iat = 0; iat < static_cast(this->Mi_.size()); iat++) + { + double occ = 0.0; + const int nh = nh_iat[iat]; + for (int ih = 0; ih < nh; ih++) + { + const int index = ib * nkb + begin_ih + ih; + occ += (conj(becp[index]) * becp[index]).real(); + } + this->Mi_[iat].z += weight * occ * sign; + begin_ih += nh; + } + } + } +} + +>>>>>>> 4c7881606 (fix(soc): reconcile DeltaSpin/tests with reverted nspin=4 m_y convention) template int SpinConstrain::get_nw() const { diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h index cdac9a615ff..37d81f26ef4 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.h +++ b/source/source_lcao/module_deltaspin/spin_constrain.h @@ -60,6 +60,36 @@ namespace spinconstrain { +<<<<<<< HEAD +======= +/** + * @brief Convert spinor occupation matrix to magnetic moment vector using Pauli matrices. + * + * @details For a two-component spinor wavefunction, the spin density matrix is: + * occ = |a|^2 conj(a)*b | = | (1+Mz)/2 (Mx+iMy)/2 | + * |conj(b)*a |b|^2 | | (Mx-iMy)/2 (1-Mz)/2 | + * NOTE occ is built conj-first (occ[1]=conj(c_up)*c_dn), i.e. the complex conjugate of the physical + * spin-density matrix, so occ[1]=(Mx+iMy)/2 and the bare Im formula below recovers the physical My. + * The magnetic moment components are extracted via Pauli matrix traces: + * Mx = Tr(rho * sigma_x) = occ[1] + occ[2] (real part) + * My = Tr(rho * sigma_y) = Im(occ[1] - occ[2]) (bare; occ is conj-first, see ref/MSG 2026-08 #7664) + * Mz = Tr(rho * sigma_z) = occ[0] - occ[3] (real part) + * where occ = {|a|^2, a*b, b*a, |b|^2} from becp coefficients. + * + * @param occ 4-element array of occupation matrix elements (complex) + * @param weight k-point weight for integration + * @return 3D magnetic moment vector (Mx, My, Mz) in Bohr magnetons + */ +inline ModuleBase::Vector3 pauli_to_moment(const std::complex occ[4], double weight) +{ + return ModuleBase::Vector3( + weight * (occ[1] + occ[2]).real(), + weight * (occ[1] - occ[2]).imag(), + weight * (occ[0] - occ[3]).real() + ); +} + +>>>>>>> 4c7881606 (fix(soc): reconcile DeltaSpin/tests with reverted nspin=4 m_y convention) struct ScAtomData; /** diff --git a/source/source_lcao/module_deltaspin/test/deltaspin_core_test.cpp b/source/source_lcao/module_deltaspin/test/deltaspin_core_test.cpp index 4cf6b4fa99d..30c036eb4e8 100644 --- a/source/source_lcao/module_deltaspin/test/deltaspin_core_test.cpp +++ b/source/source_lcao/module_deltaspin/test/deltaspin_core_test.cpp @@ -23,15 +23,17 @@ struct Vec3i { int x, y, z; }; // 1. pauli_to_moment: spinor -> magnetic moment // // Mx = w * (occ[1] + occ[2]).real() -// My = -w * (occ[1] - occ[2]).imag() (from sigma_y = [[0,-i],[i,0]]) +// My = w * (occ[1] - occ[2]).imag() (bare; occ is conj-first, occ[1]=conj(c_up)*c_dn, +// so occ[1]=(Mx+iMy)/2 and the bare Im recovers physical My) // Mz = w * (occ[0] - occ[3]).real() +// (mirrors spin_constrain.h::pauli_to_moment; #7664 flipped My, #7748 reverted the code sign) // ===================================================================== static Vec3 pauli_to_moment(const std::complex occ[4], double weight) { return { weight * (occ[1] + occ[2]).real(), - -weight * (occ[1] - occ[2]).imag(), + weight * (occ[1] - occ[2]).imag(), weight * (occ[0] - occ[3]).real() }; } @@ -81,11 +83,12 @@ TEST_F(PauliToMomentTest, GeneralCase_AllComponents) occ[2] = {0.1, -0.2}; // conj of occ[1] occ[3] = {0.4, 0.0}; auto M = pauli_to_moment(occ, 1.0); + // occ[1]=conj(c_up)*c_dn=(Mx+iMy)/2, so physical My=2*Im(occ[1])=0.4 // Mx = (0.1+0.2i + 0.1-0.2i).real = 0.2 - // My = -(0.1+0.2i - (0.1-0.2i)).imag = -(0+0.4i).imag = -0.4 + // My = (0.1+0.2i - (0.1-0.2i)).imag = (0+0.4i).imag = 0.4 // Mz = (0.6 - 0.4) = 0.2 EXPECT_NEAR(M.x, 0.2, 1e-15); - EXPECT_NEAR(M.y, -0.4, 1e-15); + EXPECT_NEAR(M.y, 0.4, 1e-15); EXPECT_NEAR(M.z, 0.2, 1e-15); } diff --git a/source/source_lcao/module_operator_lcao/dspin_lcao.cpp b/source/source_lcao/module_operator_lcao/dspin_lcao.cpp index 2d994ac37fd..be47e41ba45 100644 --- a/source/source_lcao/module_operator_lcao/dspin_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dspin_lcao.cpp @@ -56,10 +56,15 @@ inline void cal_coeff_lambda(const std::vector& current_lambda, std::vec coefficients[1] = -current_lambda[0]; } inline void cal_coeff_lambda(const std::vector& current_lambda, std::vector>& coefficients) -{// {\lambda^{I,3}, \lambda^{I,1}-i\lambda^{I,2}, \lambda^{I,1}+i\lambda^{I,2}, -\lambda^{I,3}} +{// {\lambda^{I,3}, \lambda^{I,1}+i\lambda^{I,2}, \lambda^{I,1}-i\lambda^{I,2}, -\lambda^{I,3}} + // The occupation matrix is built conj-first (occ[1]=conj(c_up)*c_dn, spin_constrain.cpp), + // and pauli_to_moment measures the physical +m_y from it (bare +Im). The lambda operator must + // drive that measured moment with the matching feedback sign, i.e. the pre-#7664 convention + // lambda_{ud}=lambda_x+i*lambda_y. #7664 flipped this together with the measurement (mirror-y); + // #7748 reverted the measurement but not this, leaving the constraint loop driving the y-mirror. coefficients[0] = std::complex(current_lambda[2], 0.0); - coefficients[1] = std::complex(current_lambda[0] , -current_lambda[1]); - coefficients[2] = std::complex(current_lambda[0] , current_lambda[1]); + coefficients[1] = std::complex(current_lambda[0] , current_lambda[1]); + coefficients[2] = std::complex(current_lambda[0] , -1 * current_lambda[1]); coefficients[3] = std::complex(-1 * current_lambda[2], 0.0); } From 30e892f2a10a5aebabf0fa7d00895ecfe06d1bbc Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 18 Aug 2026 10:10:08 -0400 Subject: [PATCH 5/6] test(soc): SpinConventionCoupling now exercises the real func_xyz_to_updown Replaces the local Pauli re-implementation with the actual elecstate::DensityMatrix_Tools::func_xyz_to_updown, fed the conj-first stored DM block (DM=conj(P)) exactly as the runtime does, and asserts (a) it recovers the physical magnetization and (b) spin_so3 (the rotation psymmg_soc applies to the grid) agrees with the SU(2) block rotation + real extraction. The two conventions can no longer drift apart silently; the test fails on the #7664 m_y flip. Links density_matrix.cpp only (no io/mocks needed). Co-Authored-By: Claude Opus 4.8 --- .../module_symmetry/test/CMakeLists.txt | 5 ++ .../test/symm_rho_soc_test.cpp | 76 +++++++++++-------- 2 files changed, 50 insertions(+), 31 deletions(-) diff --git a/source/source_cell/module_symmetry/test/CMakeLists.txt b/source/source_cell/module_symmetry/test/CMakeLists.txt index 62951d71139..4401d65abf6 100644 --- a/source/source_cell/module_symmetry/test/CMakeLists.txt +++ b/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 ) \ No newline at end of file diff --git a/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp b/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp index 9fcc9e458d9..01258305d67 100644 --- a/source/source_cell/module_symmetry/test/symm_rho_soc_test.cpp +++ b/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,35 +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 a spinor density -// block followed by the standard sigma_y=[[0,-i],[i,0]] Pauli decomposition: -// rho_0 = Re(uu+dd), rho_x = Re(ud+du), -// rho_y = -Im(ud) + Im(du), rho_z = Re(uu-dd). -// NOTE: this uses the PHYSICAL spin-density block (D_ud = mx - i*my), so the textbook -// -Im(ud)+Im(du) is correct here. This is NOT the frame the runtime uses: the actual -// stored DM is conj-first (DM=conj(P), cal_dm_psi), so func_xyz_to_updown / psymmg_soc -// consume conj(P) and use the BARE +Im(ud)-Im(du) (see #7832 / the DM round-trip test in -// source_estate/module_dm/test). This test is a spin_so3 sanity check in the physical -// frame; it does NOT validate the DM-path sign and must not be read as pinning the #7664 -// convention. TODO: exercise the real func_xyz_to_updown/psymmg_soc instead of this local -// re-implementation so the two conventions cannot drift apart silently. +// 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; -// PHYSICAL 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)}; } -// textbook Pauli extraction from the PHYSICAL block (distinct from func_xyz_to_updown, which -// reads the conj-first stored DM); 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 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(tmp, 0, step_trace, out); + return ModuleBase::Vector3(out[step_trace[1]], out[step_trace[2]], out[step_trace[3]]); } } // namespace @@ -238,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 mv(m[0], m[1], m[2]); - const ModuleBase::Vector3 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 mF = real_extract(stored_dm_from_phys(P)); + const ModuleBase::Vector3 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 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; } } } From 312393a5f2c6416ac30fe1f40297fe04360a7514 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 27 Aug 2026 07:27:18 -0400 Subject: [PATCH 6/6] rebase and solve conflicts --- .../module_deltaspin/spin_constrain.cpp | 86 ------------------- .../module_deltaspin/spin_constrain.h | 30 ------- 2 files changed, 116 deletions(-) diff --git a/source/source_lcao/module_deltaspin/spin_constrain.cpp b/source/source_lcao/module_deltaspin/spin_constrain.cpp index 003b2e422fd..fd9d1d83f08 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.cpp +++ b/source/source_lcao/module_deltaspin/spin_constrain.cpp @@ -130,92 +130,6 @@ int SpinConstrain::get_spin_sign(int ik) const return (this->pelec->klist->isk[ik] == 0) ? 1 : -1; } -<<<<<<< HEAD -======= -/** - * @brief Accumulate magnetic moments from projector coefficients (becp) for one k-point. - * - * @par Algorithm (npol=2, non-collinear): - * For each atom, compute the 2x2 occupation matrix from becp coefficients: - * occ[0] = sum_ih becp_up^*(ih) * becp_up(ih) = - * occ[1] = sum_ih becp_up^*(ih) * becp_dn(ih) = - * occ[2] = sum_ih becp_dn^*(ih) * becp_up(ih) = - * occ[3] = sum_ih becp_dn^*(ih) * becp_dn(ih) = - * where P_at = sum_{l,m} |alpha_{l,m}> - * Mz += weight * occ * spin_sign - * where spin_sign = +1 for spin-up, -1 for spin-down. - * - * @param becp Projector coefficients, layout: [ib * npol * nkb + spin * nkb + ih] - * @param nkb Total number of projectors across all atoms - * @param nbands Number of bands (occupied + unoccupied in the subspace) - * @param npol Number of spinor components (1 for collinear, 2 for non-collinear) - * @param ik K-point index (used for spin_sign lookup in collinear mode) - * @param wg_ik Band occupation weights for this k-point (from Fermi-Dirac) - * @param nh_iat Array of projector counts per atom: nh_iat[iat] = nproj for atom iat - */ -template -void SpinConstrain::accumulate_Mi_from_becp(const std::complex* becp, - int nkb, - int nbands, - int npol, - int ik, - const double* wg_ik, - const int* nh_iat) -{ - if (npol == 2) - { - for (int ib = 0; ib < nbands; ib++) - { - const double weight = wg_ik[ib]; - int begin_ih = 0; - for (int iat = 0; iat < static_cast(this->Mi_.size()); iat++) - { - std::complex occ[4] = {ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO}; - const int nh = nh_iat[iat]; - for (int ih = 0; ih < nh; ih++) - { - const int index = ib * 2 * nkb + begin_ih + ih; - occ[0] += conj(becp[index]) * becp[index]; - occ[1] += conj(becp[index]) * becp[index + nkb]; - occ[2] += conj(becp[index + nkb]) * becp[index]; - occ[3] += conj(becp[index + nkb]) * becp[index + nkb]; - } - this->Mi_[iat] += pauli_to_moment(occ, weight); - begin_ih += nh; - } - } - } - else // npol == 1 - { - const int sign = this->get_spin_sign(ik); - for (int ib = 0; ib < nbands; ib++) - { - const double weight = wg_ik[ib]; - int begin_ih = 0; - for (int iat = 0; iat < static_cast(this->Mi_.size()); iat++) - { - double occ = 0.0; - const int nh = nh_iat[iat]; - for (int ih = 0; ih < nh; ih++) - { - const int index = ib * nkb + begin_ih + ih; - occ += (conj(becp[index]) * becp[index]).real(); - } - this->Mi_[iat].z += weight * occ * sign; - begin_ih += nh; - } - } - } -} - ->>>>>>> 4c7881606 (fix(soc): reconcile DeltaSpin/tests with reverted nspin=4 m_y convention) template int SpinConstrain::get_nw() const { diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h index 37d81f26ef4..cdac9a615ff 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.h +++ b/source/source_lcao/module_deltaspin/spin_constrain.h @@ -60,36 +60,6 @@ namespace spinconstrain { -<<<<<<< HEAD -======= -/** - * @brief Convert spinor occupation matrix to magnetic moment vector using Pauli matrices. - * - * @details For a two-component spinor wavefunction, the spin density matrix is: - * occ = |a|^2 conj(a)*b | = | (1+Mz)/2 (Mx+iMy)/2 | - * |conj(b)*a |b|^2 | | (Mx-iMy)/2 (1-Mz)/2 | - * NOTE occ is built conj-first (occ[1]=conj(c_up)*c_dn), i.e. the complex conjugate of the physical - * spin-density matrix, so occ[1]=(Mx+iMy)/2 and the bare Im formula below recovers the physical My. - * The magnetic moment components are extracted via Pauli matrix traces: - * Mx = Tr(rho * sigma_x) = occ[1] + occ[2] (real part) - * My = Tr(rho * sigma_y) = Im(occ[1] - occ[2]) (bare; occ is conj-first, see ref/MSG 2026-08 #7664) - * Mz = Tr(rho * sigma_z) = occ[0] - occ[3] (real part) - * where occ = {|a|^2, a*b, b*a, |b|^2} from becp coefficients. - * - * @param occ 4-element array of occupation matrix elements (complex) - * @param weight k-point weight for integration - * @return 3D magnetic moment vector (Mx, My, Mz) in Bohr magnetons - */ -inline ModuleBase::Vector3 pauli_to_moment(const std::complex occ[4], double weight) -{ - return ModuleBase::Vector3( - weight * (occ[1] + occ[2]).real(), - weight * (occ[1] - occ[2]).imag(), - weight * (occ[0] - occ[3]).real() - ); -} - ->>>>>>> 4c7881606 (fix(soc): reconcile DeltaSpin/tests with reverted nspin=4 m_y convention) struct ScAtomData; /**