From 022948c056a8be2aa07e9b077e5ae024646913d0 Mon Sep 17 00:00:00 2001
From: dyzheng <zhengdy@aisi.ac.cn>
Date: Sat, 20 Jun 2026 22:12:48 +0800
Subject: [PATCH] =?UTF-8?q?fix:=20correct=20Fourier=20transform=20sign=20i?=
 =?UTF-8?q?n=20k=E2=86=92R=20density=20matrix=20transform?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The inverse Fourier transform from k-space to real-space density matrix
used e^{+ikR} instead of the correct e^{-ikR} as specified in
nao_lcao_force_stress_derivation.md Section 3.

This caused the stored DMR to be D_{μν}(-R) instead of the required
trace-DM D̄_{μν}(R) = D_{νμ}(-R), leading to incorrect nonlocal
pseudopotential forces for non-collinear (nspin=4) calculations with SOC.

Changes:
- density_matrix.cpp: Fix sign in cal_DMR(), cal_DMR_td(), cal_DMR_full()
  from TK(cosp, sinp) to TK(cosp, -sinp)
- dmr_complex.cpp: Fix sign in linear-response cal_DMR() specialization

Verified by Hermiticity check D(R) = D†(-R) at machine precision.
---
 .../module_dm/density_matrix.cpp              | 20 +++++++++----------
 .../module_lr/dm_trans/dmr_complex.cpp        |  6 +++---
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/source/source_estate/module_dm/density_matrix.cpp b/source/source_estate/module_dm/density_matrix.cpp
index 7414e00fd3b..e64e5d0bc7e 100644
--- a/source/source_estate/module_dm/density_matrix.cpp
+++ b/source/source_estate/module_dm/density_matrix.cpp
@@ -110,13 +110,13 @@ 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) = (1/Nk) * sum_k D(k) * exp(-i*k*R)
+                    // See Sec. 3 of nao_lcao_force_stress_derivation.md
                     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 +265,13 @@ 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) = (1/Nk) * sum_k D(k) * exp(-i*k*R)
+                    // See Sec. 3 of nao_lcao_force_stress_derivation.md
                     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 +422,13 @@ 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}
-                const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
+                // Inverse Fourier transform: D(R) = (1/Nk) * sum_k D(k) * exp(-i*k*R)
+                // See Sec. 3 of nao_lcao_force_stress_derivation.md
+                const ModuleBase::Vector3<double> dR(R_index.x, R_index.y, R_index.z);
                 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);
             }
         }
 
diff --git a/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp b/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp
index a973155c32c..42d191c1718 100644
--- a/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp
+++ b/source/source_lcao/module_lr/dm_trans/dmr_complex.cpp
@@ -43,13 +43,13 @@ namespace elecstate
                         for (int ik = 0; ik < this->_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) = (1/Nk) * sum_k D(k) * exp(-i*k*R)
+                            // See Sec. 3 of nao_lcao_force_stress_derivation.md
                             const ModuleBase::Vector3<double> dR(r_index[0], r_index[1], r_index[2]);
                             const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                             double sinp = 0.0, cosp = 0.0;
                             ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                            const std::complex<double> kphase = std::complex<double>(cosp, sinp);
+                            const std::complex<double> kphase = std::complex<double>(cosp, -sinp);
                             // set DMR element
                             std::complex<double>* tmp_DMR_pointer = tmp_matrix->get_pointer();
                             const std::complex<double>* tmp_DMK_pointer