DMK->DMR inverse Fourier sign (e^{-ikR}) and frozen density-matrix conventions for LCAO force/stress

[dyzheng]

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

1. Fourier sign consistency between cal_DMR*, dmr_complex.cpp, and
   folding_HR.
2. The closed-trace protection argument used by the force/stress operators
   (any half-R-set storage optimization must keep the pairing).
3. 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).