Fix: LCAO nspin=2/4 density-matrix chain (DMK -> DMR), the SOC nonlocal force/stress path, and the force/stress contraction tests

source/source_estate/module_dm/cal_dm_psi.cpp

@@ -24,6 +24,9 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
 
     // dm = wfc.T * wg * wfc.conj()
     // dm[is](iw1,iw2) = \sum_{ib} wfc[is](ib,iw1).T * wg(is,ib) * wfc[is](ib,iw2).conj()
+    // i.e. DMK(μ,ν) = \sum_ib wg_ib * C*_{μ,ib} * C_{ν,ib} = (C f C†)^T = D_std^T.
+    // Consumers that read DMK elements explicitly by (μ,ν) must follow this
+    // transposed convention, not treat DMK as C f C† itself.
 
     for (int ik = 0; ik < wfc.get_nk(); ++ik)
     {

source/source_estate/module_dm/density_matrix.cpp

@@ -110,13 +110,14 @@ void DensityMatrix_Tools::cal_DMR(
                 for(int ik = 0; ik < dm._nk; ++ik)
                 {
                     if(ik_in >= 0 && ik_in != ik) { continue; }
-                    // cal k_phase
-                    // if TK==std::complex<double>, kphase is e^{ikR}
+                    // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R)
+                    // k-point weights are embedded in DMK, so there is no 1/Nk
+                    // prefactor.
                     const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
                     const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                     double sinp, cosp;
                     ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                    kphase_vec[ik][iR] = TK(cosp, sinp);
+                    kphase_vec[ik][iR] = TK(cosp, -sinp);
                 }
             }
 
@@ -265,13 +266,14 @@ void DensityMatrix_Tools::cal_DMR_td(
                 for(int ik = 0; ik < dm._nk; ++ik)
                 {
                     if(ik_in >= 0 && ik_in != ik) { continue; }
-                    // cal k_phase
-                    // if TK==std::complex<double>, kphase is e^{ikR}
+                    // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R)
+                    // k-point weights are embedded in DMK, so there is no 1/Nk
+                    // prefactor.
                     const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
                     const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                     double sinp, cosp;
                     ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                    kphase_vec[ik][iR] = TK(cosp, sinp);
+                    kphase_vec[ik][iR] = TK(cosp, -sinp);
                     if(PARAM.inp.td_stype==2)
                     {
                         //phase for hybrid gauge tddft
@@ -422,13 +424,14 @@ void DensityMatrix_Tools::cal_DMR_full(
             for(int ik = 0; ik < dm._nk; ++ik)
             {
                 if(ik_in >= 0 && ik_in != ik) { continue; }
-                // cal k_phase
-                // if TK==std::complex<double>, kphase is e^{ikR}
+                // Inverse Fourier transform: D(R) = sum_k D(k) * exp(-i*k*R)
+                // Phase factor: exp(-i*k*R) = cos(k·R) - i*sin(k·R)
+                // k-point weights are embedded in DMK, so there is no 1/Nk prefactor.
                 const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
                 const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                 double sinp, cosp;
                 ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                kphase_vec[ik][iR] = TK(cosp, sinp);
+                kphase_vec[ik][iR] = TK(cosp, -sinp);
             }
         }
 

source/source_estate/module_dm/density_matrix.h

@@ -9,6 +9,26 @@
 
 namespace elecstate
 {
+// ---------------------------------------------------------------------------
+// Density-matrix conventions (DMK/DMR), frozen by the DM/DMK/DMR chain.
+//
+// 1. DMK(μ,ν;k) = Σ_n f_nk C*_{μn}(k) C_{νn}(k) = (C f C†)^T = D_std^T.
+//    μ = row (bra) orbital, ν = column (ket) orbital. Stored as the transpose
+//    of the textbook D_std (cal_dm_psi.cpp); for Hermitian D this equals the
+//    conjugate. Consumers that read DMK elements explicitly by (μ,ν) must
+//    follow this convention.
+// 2. D(R) = Σ_k e^{-ikR} DMK(k), k·R in direct coordinates with 2π.
+//    k-point weights w_k are embedded in DMK (wg = w_k*occ), so cal_DMR must
+//    NOT multiply by 1/Nk again. The forward pair used by operators is
+//    O(k) = Σ_R e^{+ikR} O(R) (folding_HR).
+// 3. Closed-trace protection: a contraction Σ_{μνR} D(μ,ν;R)·O(μ,ν;R) equals
+//    Σ_k w_k·Re Tr(DMK(k)·O(k)) iff
+//    (a) the contraction is a closed Frobenius trace over all orbitals,
+//    (b) the paired operator is Hermitian per k (the HContainer stores
+//        (iat1,iat2,R) and (iat2,iat1,-R) with O(-R) = O(R)†),
+//    (c) only the real part is kept.
+// ---------------------------------------------------------------------------
+
 /**
  * @brief DensityMatrix Class
  * <TK,TR> = <double,double> for Gamma-only calculation

source/source_estate/module_dm/test/CMakeLists.txt

@@ -46,5 +46,6 @@ AddTest(
   ${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_hamilt/module_hcontainer/func_folding.cpp
   ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp
 )