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