Describe the bug
1. Problem
The real-space density matrix DMR is produced from the k-space density
matrix DMK by an inverse Fourier transform. In ABACUS the Hamiltonian is
folded with e^{+ikR} (folding_HR, func_folding.cpp); consistency of the
SCF loop and of any force/stress contraction therefore requires the inverse
transform to use the opposite phase, e^{-ikR}.
The LCAO code did not consistently satisfy this pairing:
- cal_DMR / cal_DMR_td / cal_DMR_full in
source/source_estate/module_dm/density_matrix.cpp and the linear-response
path dmr_complex.cpp used a k-phase whose sign disagreed with folding_HR,
so the (iat1, iat2, R) pairing of DMR did not match the R-pairing used by
the nonlocal force/stress operators. The mismatch corrupts the nspin=4
(SOC) nonlocal force, where the same R vector enters both the overlap and
the derivative of the beta projectors.
- The underlying density-matrix conventions were undocumented and unfrozen:
DMK is stored transposed (DMK = (C f C^+)^T, built by cal_dm_psi with
the conjugate applied to the first index), k-point weights are embedded in
DMK (no explicit 1/Nk), and the force/stress trace is protected by the
closed trace over all R directions. Without a written convention, any change
to the Fourier phase or to the storage layout silently breaks forces.
- operator_fs_utils contained stale "factor of 2" comments that no longer
matched the full-R-set iteration, and the Mulliken M_y inline comment did
not match its algebra (code unchanged).
2. Proposed changes
2.1 Fourier sign fix (source)
- cal_DMR, cal_DMR_td, cal_DMR_full (density_matrix.cpp): use
e^{-ikR} as the inverse transform of folding_HR's e^{+ikR}.
- Linear response dmr_complex.cpp and exciton_plotter.cpp: same sign
convention for the LR DM path.
- density_matrix.h: document the frozen conventions in one block:
- DMK(mu,nu;k) = sum_n f_nk C*_{mun}(k) C_{nun}(k) = (C f C^+)^T = D_std^T;
- k weights embedded in DMK, no 1/Nk factor;
- DMR(iat2,iat1,-R) pairing and the closed-trace protection condition.
- operator_fs_utils: replace the stale factor-of-2 comments; the full-R-set
iteration uses factor 1.0.
- output_mulliken.cpp: fix the M_y inline-comment algebra (code unchanged).
2.2 Unit tests (source/source_estate/module_dm/test)
- T1 Fourier round-trip: cal_DMR (e^{-ikR}) vs folding_HR (e^{+ikR}),
with a sign sentinel (a non-Gamma, non-TRIM k grid). Flipping the phase back
to +sinp turns the test red.
- T2 DMR Hermiticity: D(iat2,iat1,-R) equals D(iat1,iat2,R)^T (real) or
^dagger (complex, nspin=4).
- T3 closed-trace equivalence, including a non-symmetric-in-R operator.
- T8 full-direction pairing storage guard: every (iat1,iat2,R) block must
have a (iat2,iat1,-R) counterpart (the force/stress paths rely on it).
- test_dm_trace.cpp (source/source_lcao/module_operator_lcao/test):
closed-trace W_k = sum_k Re Tr( DMK(k) O(k) ) checks for the
overlap/kinetic derivative operators.
2.3 Integration tests
- T4 240_NO_KP_15_SO_FD (nspin=4, multi-k, SOC, GaAs): check_extra.py
compares the analytic nonlocal force against a central finite difference
(delta = 1e-3 Bohr, threshold 1e-4 eV/Bohr) and checks the acoustic sum
rule. Measured max deviation ~5.8e-5 eV/Bohr.
- T4 260_NO_DJ_PK_PU_SO_4K (nspin=4, DFT+U, 4x1x1 grid with non-TRIM k
points): reference regenerated for e^{-ikR} (etot -6791.4650419973 eV).
Flipping the sign back to +sinp shifts etot by ~0.44 eV and the Fe1 force
by ~4 eV/Ang, so the case turns red under the wrong sign.
- T5 18_NO_hyb_TDDFT (nspin=2, hybrid TDDFT): check_extra.py physics
checks (zero-field current ~5e-17, full-k vs per-k DMR path deviation
~3e-18, +/-k cancellation = 0).
- T6 260_NO_DJ_PK_PU_SO_4K_MUL (nspin=4 multi-k Mulliken): Mulliken
charge/magnetization cross-checked against the rho-cell magnetization
(max |dev| 0.044 uB < 0.1 uB). No anomalous m_y deviation present.
- T7 260_NO_DJ_PK_PU_AFM_URAMPING_MUL (nspin=2 AFM): opposite-sign Fe
moments +1.92/-3.73 uB, spin1+spin2 = total charge 44, small O moments,
Mulliken sum vs rho-cell magnetization.
- Autotest.sh: check_out now uses continue instead of break after
totaltimeref (so case-local check_extra.py keys are actually compared)
and runs check_extra.py right after catch_properties.sh, appending its
output to result.out/result.ref.
3. Verification
- OMP_NUM_THREADS=1, mpirun -np 1: full CASES_CPU.txt run passes 27/27
keys; T5 passes 6/6. All new cases converge (SCF CONVERGED, guarded by
*_converged_pass keys).
- Unit suite: ctest -R 'ESTATE_dm|LCAO_operator' passes (T1-T8); the
parallel_operator_tests.sh failure is a pre-existing build-dir artifact.
- Sign sensitivity: flipping -sinp back to +sinp turns T1 red and shifts
the T4 DFT+U reference (see 2.3).
4. Notes for the PR
- The derivation document that accompanied the fix
(docs/nao_lcao_force_stress_derivation.md) has been dropped from the
branch; the conventions now live in the density_matrix.h comment block and
the T1/T2/T3/T8 tests.
- The change is independent of the nspin=4 m_y convention discussion
(PR Fix: revert #7664 $m_y$ convention #7832): this branch locks the real-space Fourier phase; Fix: revert #7664 $m_y$ convention #7832 concerns
the spin-space rho_y sign. The two patches are orthogonal (no overlapping
hunks in density_matrix.cpp).
5. Requested review focus
- Fourier sign consistency between cal_DMR*, dmr_complex.cpp, and
folding_HR.
- The closed-trace protection argument used by the force/stress operators
(any half-R-set storage optimization must keep the pairing).
- Whether the Mulliken M_y comment correction should be accompanied by a
code fix in a follow-up (currently comment-only, code unchanged).
Expected behavior
..
To Reproduce
No response
Environment
No response
Additional Context
No response
Task list for Issue attackers (only for developers)
- Verify the issue is not a duplicate.
- Describe the bug.
- Steps to reproduce.
- Expected behavior.
- Error message.
- Environment details.
- Additional context.
- Assign a priority level (low, medium, high, urgent).
- Assign the issue to a team member.
- Label the issue with relevant tags.
- Identify possible related issues.
- Create a unit test or automated test to reproduce the bug (if applicable).
- Fix the bug.
- Test the fix.
- Update documentation (if necessary).
- Close the issue and inform the reporter (if applicable).
Reactions are currently unavailable
Describe the bug
1. Problem
The real-space density matrix DMR is produced from the k-space density
matrix DMK by an inverse Fourier transform. In ABACUS the Hamiltonian is
folded with e^{+ikR} (folding_HR, func_folding.cpp); consistency of the
SCF loop and of any force/stress contraction therefore requires the inverse
transform to use the opposite phase, e^{-ikR}.
The LCAO code did not consistently satisfy this pairing:
source/source_estate/module_dm/density_matrix.cpp and the linear-response
path dmr_complex.cpp used a k-phase whose sign disagreed with folding_HR,
so the (iat1, iat2, R) pairing of DMR did not match the R-pairing used by
the nonlocal force/stress operators. The mismatch corrupts the nspin=4
(SOC) nonlocal force, where the same R vector enters both the overlap and
the derivative of the beta projectors.
DMK is stored transposed (DMK = (C f C^+)^T, built by cal_dm_psi with
the conjugate applied to the first index), k-point weights are embedded in
DMK (no explicit 1/Nk), and the force/stress trace is protected by the
closed trace over all R directions. Without a written convention, any change
to the Fourier phase or to the storage layout silently breaks forces.
matched the full-R-set iteration, and the Mulliken M_y inline comment did
not match its algebra (code unchanged).
2. Proposed changes
2.1 Fourier sign fix (source)
e^{-ikR} as the inverse transform of folding_HR's e^{+ikR}.
convention for the LR DM path.
iteration uses factor 1.0.
2.2 Unit tests (source/source_estate/module_dm/test)
with a sign sentinel (a non-Gamma, non-TRIM k grid). Flipping the phase back
to +sinp turns the test red.
^dagger (complex, nspin=4).
have a (iat2,iat1,-R) counterpart (the force/stress paths rely on it).
closed-trace W_k = sum_k Re Tr( DMK(k) O(k) ) checks for the
overlap/kinetic derivative operators.
2.3 Integration tests
compares the analytic nonlocal force against a central finite difference
(delta = 1e-3 Bohr, threshold 1e-4 eV/Bohr) and checks the acoustic sum
rule. Measured max deviation ~5.8e-5 eV/Bohr.
points): reference regenerated for e^{-ikR} (etot -6791.4650419973 eV).
Flipping the sign back to +sinp shifts etot by ~0.44 eV and the Fe1 force
by ~4 eV/Ang, so the case turns red under the wrong sign.
checks (zero-field current ~5e-17, full-k vs per-k DMR path deviation
~3e-18, +/-k cancellation = 0).
charge/magnetization cross-checked against the rho-cell magnetization
(max |dev| 0.044 uB < 0.1 uB). No anomalous m_y deviation present.
moments +1.92/-3.73 uB, spin1+spin2 = total charge 44, small O moments,
Mulliken sum vs rho-cell magnetization.
totaltimeref (so case-local check_extra.py keys are actually compared)
and runs check_extra.py right after catch_properties.sh, appending its
output to result.out/result.ref.
3. Verification
keys; T5 passes 6/6. All new cases converge (SCF CONVERGED, guarded by
*_converged_pass keys).
parallel_operator_tests.sh failure is a pre-existing build-dir artifact.
the T4 DFT+U reference (see 2.3).
4. Notes for the PR
(docs/nao_lcao_force_stress_derivation.md) has been dropped from the
branch; the conventions now live in the density_matrix.h comment block and
the T1/T2/T3/T8 tests.
(PR Fix: revert #7664 $m_y$ convention #7832): this branch locks the real-space Fourier phase; Fix: revert #7664 $m_y$ convention #7832 concerns
the spin-space rho_y sign. The two patches are orthogonal (no overlapping
hunks in density_matrix.cpp).
5. Requested review focus
folding_HR.
(any half-R-set storage optimization must keep the pairing).
code fix in a follow-up (currently comment-only, code unchanged).
Expected behavior
..
To Reproduce
No response
Environment
No response
Additional Context
No response
Task list for Issue attackers (only for developers)