From d5358d5b40cf45ed5525190a2629ff4ff3505cb0 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Sun, 13 Oct 2024 14:48:29 +0800
Subject: [PATCH 1/6] remove dependence of DensityMatrix on K_Vectors

---
 source/module_elecstate/elecstate_lcao.cpp    |  3 +-
 .../module_dm/density_matrix.cpp              | 73 +++++--------------
 .../module_dm/density_matrix.h                | 19 +++--
 .../module_dm/test/test_cal_dm_R.cpp          |  8 +-
 .../module_dm/test/test_dm_R_init.cpp         | 16 ++--
 .../module_dm/test/test_dm_constructor.cpp    |  4 +-
 .../module_dm/test/test_dm_io.cpp             |  4 +-
 .../hamilt_lcaodft/hamilt_lcao.h              |  1 +
 source/module_io/get_pchg_lcao.cpp            |  3 +-
 source/module_io/get_pchg_lcao.h              |  1 +
 source/module_io/td_current_io.cpp            |  7 +-
 source/module_lr/dm_trans/dmr_complex.cpp     |  2 +-
 source/module_lr/hamilt_casida.h              |  2 +-
 source/module_lr/lr_spectrum.cpp              |  6 +-
 source/module_lr/lr_spectrum.h                |  1 +
 .../operator_casida/operator_lr_hxc.cpp       |  2 +-
 .../operator_casida/operator_lr_hxc.h         |  3 +-
 17 files changed, 63 insertions(+), 92 deletions(-)

diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp
index 0d11153183..fdcbac0e43 100644
--- a/source/module_elecstate/elecstate_lcao.cpp
+++ b/source/module_elecstate/elecstate_lcao.cpp
@@ -131,7 +131,8 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
 template <typename TK>
 void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
 {
-    this->DM = new DensityMatrix<TK, double>(kv, paraV, nspin);
+    const int nspin_dm = std::map<int, int>({ {1, 1}, {2, 2}, {4, 1} })[nspin];
+    this->DM = new DensityMatrix<TK, double>(paraV, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
 }
 
 template <>
diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp
index 22152bb093..8ff877a741 100644
--- a/source/module_elecstate/module_dm/density_matrix.cpp
+++ b/source/module_elecstate/module_dm/density_matrix.cpp
@@ -24,62 +24,24 @@ DensityMatrix<TK, TR>::~DensityMatrix()
     }
 }
 
-// constructor for multi-k
 template <typename TK, typename TR>
-DensityMatrix<TK, TR>::DensityMatrix(const K_Vectors* kv_in, const Parallel_Orbitals* paraV_in, const int nspin)
+DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin, const std::vector<ModuleBase::Vector3<double>>& kvec_d, const int nk)
+    : _paraV(paraV_in), _nspin(nspin), _kvec_d(kvec_d), _nks((nk > 0 && nk <= _kvec_d.size()) ? nk : _kvec_d.size())
 {
     ModuleBase::TITLE("DensityMatrix", "DensityMatrix-MK");
-    this->_kv = kv_in;
-    this->_paraV = paraV_in;
-    // set this->_nspin
-    if (nspin == 1 || nspin == 4)
-    {
-        this->_nspin = 1;
-    }
-    else if (nspin == 2)
-    {
-        this->_nspin = 2;
-#ifdef __DEBUG
-        assert(kv_in->get_nks() % 2 == 0);
-#endif
-    }
-    else
-    {
-        throw std::string("nspin must be 1, 2 or 4");
-    }
-    // set this->_nks, which is real number of k-points
-    this->_nks = kv_in->get_nks() / this->_nspin;
-    // allocate memory for _DMK
-    this->_DMK.resize(this->_kv->get_nks());
-    for (int ik = 0; ik < this->_kv->get_nks(); ik++)
+    const int nks = _nks * _nspin;
+    this->_DMK.resize(nks);
+    for (int ik = 0; ik < nks; ik++)
     {
         this->_DMK[ik].resize(this->_paraV->get_row_size() * this->_paraV->get_col_size());
     }
     ModuleBase::Memory::record("DensityMatrix::DMK", this->_DMK.size() * this->_DMK[0].size() * sizeof(TK));
 }
 
-// constructor for Gamma-Only
 template <typename TK, typename TR>
-DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin)
+DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin) :_paraV(paraV_in), _nspin(nspin), _kvec_d({ ModuleBase::Vector3<double>(0,0,0) }), _nks(1)
 {
     ModuleBase::TITLE("DensityMatrix", "DensityMatrix-GO");
-    this->_paraV = paraV_in;
-    // set this->_nspin
-    if (nspin == 1 || nspin == 4)
-    {
-        this->_nspin = 1;
-    }
-    else if (nspin == 2)
-    {
-        this->_nspin = 2;
-    }
-    else
-    {
-        throw std::string("nspin must be 1, 2 or 4");
-    }
-    // set this->_nks, which is real number of k-points
-    this->_nks = 1;
-    // allocate memory for _DMK
     this->_DMK.resize(_nspin);
     for (int ik = 0; ik < this->_nspin; ik++)
     {
@@ -328,10 +290,9 @@ template <typename TK, typename TR>
 int DensityMatrix<TK, TR>::get_DMK_nks() const
 {
 #ifdef __DEBUG
-    assert(this->_DMK.size() != 0);
-    assert(this->_kv != nullptr);
+    assert(this->_DMK.size() == _nks * _nspin);
 #endif
-    return this->_kv->get_nks();
+    return _nks * _nspin;
 }
 
 template <typename TK, typename TR>
@@ -440,7 +401,7 @@ void DensityMatrix<TK, TR>::cal_DMR_test()
                     // 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]);
-                    const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                    const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                     double sinp, cosp;
                     ModuleBase::libm::sincos(arg, &sinp, &cosp);
                     std::complex<double> kphase = std::complex<double>(cosp, sinp);
@@ -515,7 +476,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
                         // 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]);
-                        const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                        const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                         double sinp, cosp;
                         ModuleBase::libm::sincos(arg, &sinp, &cosp);
                         std::complex<double> kphase = std::complex<double>(cosp, sinp);
@@ -561,7 +522,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
                         // 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]);
-                        const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                        const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                         double sinp, cosp;
                         ModuleBase::libm::sincos(arg, &sinp, &cosp);
                         std::complex<double> kphase = std::complex<double>(cosp, sinp);
@@ -678,7 +639,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik)
                     // 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]);
-                    const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                    const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                     double sinp, cosp;
                     ModuleBase::libm::sincos(arg, &sinp, &cosp);
                     std::complex<double> kphase = std::complex<double>(cosp, sinp);
@@ -855,9 +816,9 @@ void DensityMatrix<TK, TR>::read_DMK(const std::string directory, const int ispi
         // quit the program or not.
         bool quit = false;
 
-        ModuleBase::CHECK_DOUBLE(ifs, this->_kv->kvec_d[ik].x, quit);
-        ModuleBase::CHECK_DOUBLE(ifs, this->_kv->kvec_d[ik].y, quit);
-        ModuleBase::CHECK_DOUBLE(ifs, this->_kv->kvec_d[ik].z, quit);
+        ModuleBase::CHECK_DOUBLE(ifs, this->_kvec_d[ik].x, quit);
+        ModuleBase::CHECK_DOUBLE(ifs, this->_kvec_d[ik].y, quit);
+        ModuleBase::CHECK_DOUBLE(ifs, this->_kvec_d[ik].z, quit);
         ModuleBase::CHECK_INT(ifs, this->_paraV->nrow);
         ModuleBase::CHECK_INT(ifs, this->_paraV->ncol);
     } // If file exist, read in data.
@@ -890,7 +851,7 @@ void DensityMatrix<double, double>::write_DMK(const std::string directory, const
     {
         ModuleBase::WARNING("elecstate::write_dmk", "Can't create DENSITY MATRIX File!");
     }
-    ofs << this->_kv->kvec_d[ik].x << " " << this->_kv->kvec_d[ik].y << " " << this->_kv->kvec_d[ik].z << std::endl;
+    ofs << this->_kvec_d[ik].x << " " << this->_kvec_d[ik].y << " " << this->_kvec_d[ik].z << std::endl;
     ofs << "\n  " << this->_paraV->nrow << " " << this->_paraV->ncol << std::endl;
 
     ofs << std::setprecision(3);
@@ -927,7 +888,7 @@ void DensityMatrix<std::complex<double>, double>::write_DMK(const std::string di
     {
         ModuleBase::WARNING("elecstate::write_dmk", "Can't create DENSITY MATRIX File!");
     }
-    ofs << this->_kv->kvec_d[ik].x << " " << this->_kv->kvec_d[ik].y << " " << this->_kv->kvec_d[ik].z << std::endl;
+    ofs << this->_kvec_d[ik].x << " " << this->_kvec_d[ik].y << " " << this->_kvec_d[ik].z << std::endl;
     ofs << "\n  " << this->_paraV->nrow << " " << this->_paraV->ncol << std::endl;
 
     ofs << std::setprecision(3);
diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h
index 0bc28c6cad..02f389099a 100644
--- a/source/module_elecstate/module_dm/density_matrix.h
+++ b/source/module_elecstate/module_dm/density_matrix.h
@@ -3,7 +3,6 @@
 
 #include <string>
 
-#include "module_cell/klist.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h"
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
@@ -45,16 +44,20 @@ class DensityMatrix
 
     /**
      * @brief Constructor of class DensityMatrix for multi-k calculation
-     * @param _kv pointer of K_Vectors object
      * @param _paraV pointer of Parallel_Orbitals object
-     * @param nspin spin setting (1 - none spin; 2 - spin; 4 - SOC)
+     * @param nspin number of spin of the density matrix, set by user according to global nspin
+     *  (usually {nspin_global -> nspin_dm} = {1->1, 2->2, 4->1}, but sometimes 2->1 like in LR-TDDFT)
+     * @param kvec_d direct coordinates of kpoints
+     * @param nk number of k-points, not always equal to K_Vectors::get_nks()/nspin_dm.
+     * if remains default or large than kvec_d.size(), it will be set to kvec_d.size()
      */
-    DensityMatrix(const K_Vectors* _kv, const Parallel_Orbitals* _paraV, const int nspin);
+    DensityMatrix(const Parallel_Orbitals* _paraV, const int nspin, const std::vector<ModuleBase::Vector3<double>>& kvec_d, const int nk = -1);
 
     /**
      * @brief Constructor of class DensityMatrix for gamma-only calculation, where kvector is not required
      * @param _paraV pointer of Parallel_Orbitals object
-     * @param nspin spin setting (1 - none spin; 2 - spin; 4 - SOC)
+     * @param nspin number of spin of the density matrix, set by user according to global nspin
+     *  (usually {nspin_global -> nspin_dm} = {1->1, 2->2, 4->1}, but sometimes 2->1 like in LR-TDDFT)
      */
     DensityMatrix(const Parallel_Orbitals* _paraV, const int nspin);
 
@@ -169,7 +172,7 @@ class DensityMatrix
      */
     const Parallel_Orbitals* get_paraV_pointer() const {return this->_paraV;}
 
-    const K_Vectors* get_kv_pointer() const {return this->_kv;}
+    const std::vector<ModuleBase::Vector3<double>>& get_kvec_d() const { return this->_kvec_d; }
 
     /**
      * @brief calculate density matrix DMR from dm(k) using blas::axpy
@@ -249,7 +252,7 @@ class DensityMatrix
     /**
      * @brief K_Vectors object, which is used to get k-point information
      */
-    const K_Vectors* _kv;
+    const std::vector<ModuleBase::Vector3<double>> _kvec_d;
 
     /**
      * @brief Parallel_Orbitals object, which contain all information of 2D block cyclic distribution
@@ -268,7 +271,7 @@ class DensityMatrix
      * _nks is not equal to _kv->get_nks() when spin-polarization is considered
      * _nks = kv->_nks / nspin
      */
-    int _nks = 0;
+    int _nks = 0;   // comment by lunasea: I want to rename it as nk
 
 
 };
diff --git a/source/module_elecstate/module_dm/test/test_cal_dm_R.cpp b/source/module_elecstate/module_dm/test/test_cal_dm_R.cpp
index 8cb0ab7539..c0802a4942 100644
--- a/source/module_elecstate/module_dm/test/test_cal_dm_R.cpp
+++ b/source/module_elecstate/module_dm/test/test_cal_dm_R.cpp
@@ -74,7 +74,7 @@ class DMTest : public testing::Test
         init_parav();
     }
 
-    void TearDown()
+    void TearDown() override
     {
         delete paraV;
         delete[] ucell.atoms[0].tau;
@@ -129,7 +129,7 @@ TEST_F(DMTest, cal_DMR_test)
     kv->set_nks(nks);
     kv->kvec_d.resize(nks);
     // construct DM
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     // set this->_DMK
     for (int is = 1; is <= nspin; is++)
     {
@@ -198,7 +198,7 @@ TEST_F(DMTest, cal_DMR_blas_double)
     kv->set_nks(nks);
     kv->kvec_d.resize(nks);
     // construct DM
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     // set this->_DMK
     for (int is = 1; is <= nspin; is++)
     {
@@ -269,7 +269,7 @@ TEST_F(DMTest, cal_DMR_blas_complex)
     kv->kvec_d[1].x = 0.5;
     kv->kvec_d[3].x = 0.5;
     // construct DM
-    elecstate::DensityMatrix<std::complex<double>, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     // set this->_DMK
     for (int is = 1; is <= nspin; is++)
     {
diff --git a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp
index 72aabf01ec..e32308d456 100644
--- a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp
+++ b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp
@@ -74,7 +74,7 @@ class DMTest : public testing::Test
         init_parav();
     }
 
-    void TearDown()
+    void TearDown() override
     {
         delete paraV;
         delete[] ucell.atoms[0].tau;
@@ -118,7 +118,7 @@ TEST_F(DMTest, DMInit1)
     // construct DM
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
     // initialize this->_DMR
     Grid_Driver gd(0,0);
     DM.init_DMR(&gd, &ucell);
@@ -145,7 +145,7 @@ TEST_F(DMTest, DMInit2)
     // construct DM
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
     // initialize Record_adj using Grid_Driver
     Grid_Driver gd(0,0);
     Record_adj ra;
@@ -208,12 +208,12 @@ TEST_F(DMTest, DMInit3)
     kv->kvec_d[1].x = 0.5;
     kv->kvec_d[3].x = 0.5;
     // construct a DM
-    elecstate::DensityMatrix<std::complex<double>, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     Grid_Driver gd(0, 0);
     DM.init_DMR(&gd, &ucell);
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     // construct another DM
-    elecstate::DensityMatrix<std::complex<double>, double> DM1(kv, paraV, nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM1(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     DM1.init_DMR(*DM.get_DMR_pointer(1));
     // compare
     EXPECT_EQ(DM1.get_DMR_pointer(2)->size_atom_pairs(), test_size * test_size);
@@ -266,7 +266,7 @@ TEST_F(DMTest, DMInit4)
         }
     }
     // construct a DM from this HContainer
-    elecstate::DensityMatrix<std::complex<double>, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     DM.init_DMR(*tmp_DMR);
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     // compare
@@ -292,11 +292,11 @@ TEST_F(DMTest, saveDMR)
     kv->kvec_d[1].x = 0.5;
     kv->kvec_d[3].x = 0.5;
     // construct a DM
-    elecstate::DensityMatrix<std::complex<double>, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     Grid_Driver gd(0, 0);
     DM.init_DMR(&gd, &ucell);
     // construct another DM
-    elecstate::DensityMatrix<std::complex<double>, double> DM_test(kv, paraV, nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM_test(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     DM_test.init_DMR(*DM.get_DMR_pointer(1));
     DM_test.save_DMR();
     EXPECT_EQ(DM_test.get_DMR_pointer(1)->get_nnr(), DM.get_DMR_pointer(1)->get_nnr());
diff --git a/source/module_elecstate/module_dm/test/test_dm_constructor.cpp b/source/module_elecstate/module_dm/test/test_dm_constructor.cpp
index f055c8b8de..c2c9832c44 100644
--- a/source/module_elecstate/module_dm/test/test_dm_constructor.cpp
+++ b/source/module_elecstate/module_dm/test/test_dm_constructor.cpp
@@ -127,7 +127,7 @@ TEST_F(DMTest, DMConstructor_nspin1)
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
     int nspin = 1;
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
     // compare
     EXPECT_EQ(DM.get_DMK_nks(), kv->get_nks());
     EXPECT_EQ(DM.get_DMK_nrow(), paraV->nrow);
@@ -196,7 +196,7 @@ TEST_F(DMTest, DMConstructor_nspin2)
     // construct DM
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks() / nspin);
     // compare
     EXPECT_EQ(DM.get_DMK_nks(), kv->get_nks());
     EXPECT_EQ(DM.get_DMK_nrow(), paraV->nrow);
diff --git a/source/module_elecstate/module_dm/test/test_dm_io.cpp b/source/module_elecstate/module_dm/test/test_dm_io.cpp
index 87bf29b937..8622bbc649 100644
--- a/source/module_elecstate/module_dm/test/test_dm_io.cpp
+++ b/source/module_elecstate/module_dm/test/test_dm_io.cpp
@@ -142,7 +142,7 @@ TEST_F(DMTest, DMConstructor1)
     int nspin = 1;
     // construct DM
     std::cout << paraV->nrow << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
     // read DMK
     std::string directory = "./support/";
     for (int is = 1; is <= nspin; ++is)
@@ -162,7 +162,7 @@ TEST_F(DMTest, DMConstructor1)
         }
     }
     // construct a new DM
-    elecstate::DensityMatrix<double, double> DM1(kv, paraV, nspin);
+    elecstate::DensityMatrix<double, double> DM1(paraV, nspin, kv->kvec_d);
     directory = "./support/output";
     for (int is = 1; is <= nspin; ++is)
     {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h
index 740fff3008..c11d419c3d 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h
@@ -2,6 +2,7 @@
 #define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_HAMILT_LCAO_H
 
 #include "module_basis/module_nao/two_center_bundle.h"
+#include "module_cell/klist.h"
 #include "module_elecstate/module_dm/density_matrix.h"
 #include "module_elecstate/potentials/potential_new.h"
 #include "module_hamilt_general/hamilt.h"
diff --git a/source/module_io/get_pchg_lcao.cpp b/source/module_io/get_pchg_lcao.cpp
index 2858147353..7910f22584 100644
--- a/source/module_io/get_pchg_lcao.cpp
+++ b/source/module_io/get_pchg_lcao.cpp
@@ -219,7 +219,8 @@ void IState_Charge::begin(Gint_k& gk,
         if (bands_picked_[ib])
         {
             // Using new density matrix inplementation (multi-k)
-            elecstate::DensityMatrix<std::complex<double>, double> DM(&kv, this->ParaV, nspin);
+            const int nspin_dm = std::map<int, int>({ {1,1},{2,2},{4,1} })[nspin];
+            elecstate::DensityMatrix<std::complex<double>, double> DM(this->ParaV, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm);
 
 #ifdef __MPI
             this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv, if_separate_k);
diff --git a/source/module_io/get_pchg_lcao.h b/source/module_io/get_pchg_lcao.h
index 57617280c8..a923397c32 100644
--- a/source/module_io/get_pchg_lcao.h
+++ b/source/module_io/get_pchg_lcao.h
@@ -1,6 +1,7 @@
 #ifndef ISTATE_CHARGE_H
 #define ISTATE_CHARGE_H
 #include "module_basis/module_pw/pw_basis.h"
+#include "module_cell/klist.h"
 #include "module_elecstate/module_dm/density_matrix.h"
 #include "module_hamilt_lcao/module_gint/gint.h"
 #include "module_hamilt_lcao/module_gint/gint_gamma.h"
diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp
index 72cd9fd412..50fae430c2 100644
--- a/source/module_io/td_current_io.cpp
+++ b/source/module_io/td_current_io.cpp
@@ -64,7 +64,7 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix<std::complex<double>, double>
                 // cal k_phase
                 // if TK==std::complex<double>, kphase is e^{ikR}
                 const ModuleBase::Vector3<double> dR(r_index.x, r_index.y, r_index.z);
-                const double arg = (DM_real.get_kv_pointer()->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                const double arg = (DM_real.get_kvec_d()[ik] * dR) * ModuleBase::TWO_PI;
                 double sinp, cosp;
                 ModuleBase::libm::sincos(arg, &sinp, &cosp);
                 std::complex<double> kphase = std::complex<double>(cosp, sinp);
@@ -156,8 +156,9 @@ void ModuleIO::write_current(const int istep,
     // construct a DensityMatrix object
     // Since the function cal_dm_psi do not suport DMR in complex type, I replace it with two DMR in double type. Should
     // be refactored in the future.
-    elecstate::DensityMatrix<std::complex<double>, double> DM_real(&kv, pv, PARAM.inp.nspin);
-    elecstate::DensityMatrix<std::complex<double>, double> DM_imag(&kv, pv, PARAM.inp.nspin);
+    const int nspin_dm = std::map<int, int>({ {1,1},{2,2},{4,1} })[PARAM.inp.nspin];
+    elecstate::DensityMatrix<std::complex<double>, double> DM_real(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm);
+    elecstate::DensityMatrix<std::complex<double>, double> DM_imag(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm);
     // calculate DMK
     elecstate::cal_dm_psi(DM_real.get_paraV_pointer(), pelec->wg, psi[0], DM_real);
 
diff --git a/source/module_lr/dm_trans/dmr_complex.cpp b/source/module_lr/dm_trans/dmr_complex.cpp
index cdc3e5139c..265e962546 100644
--- a/source/module_lr/dm_trans/dmr_complex.cpp
+++ b/source/module_lr/dm_trans/dmr_complex.cpp
@@ -48,7 +48,7 @@ namespace elecstate
                             // 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]);
-                            const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                            const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                             double sinp, cosp;
                             ModuleBase::libm::sincos(arg, &sinp, &cosp);
                             std::complex<double> kphase = std::complex<double>(cosp, sinp);
diff --git a/source/module_lr/hamilt_casida.h b/source/module_lr/hamilt_casida.h
index f203316dff..9d9e0d3ebf 100644
--- a/source/module_lr/hamilt_casida.h
+++ b/source/module_lr/hamilt_casida.h
@@ -44,7 +44,7 @@ namespace LR
             if (ri_hartree_benchmark != "aims") { assert(aims_nbasis.empty()); }
             this->classname = "HamiltCasidaLR";
             this->DM_trans.resize(1);
-            this->DM_trans[0] = LR_Util::make_unique<elecstate::DensityMatrix<T, T>>(&kv_in, pmat_in, nspin);
+            this->DM_trans[0] = LR_Util::make_unique<elecstate::DensityMatrix<T, T>>(pmat_in, nspin, kv_in.kvec_d, nk);
             // add the diag operator  (the first one)
             this->ops = new OperatorLRDiag<T>(eig_ks, pX_in, nk, nocc, nvirt);
             //add Hxc operator
diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp
index f39940bc2c..dcfefee4ea 100644
--- a/source/module_lr/lr_spectrum.cpp
+++ b/source/module_lr/lr_spectrum.cpp
@@ -33,7 +33,7 @@ void LR::LR_Spectrum<double>::oscillator_strength()
     osc.resize(X.get_nbands(), 0.0);
     // const int nspin0 = (this->nspin == 2) ? 2 : 1;   use this in NSPIN=4 implementation
     double osc_tot = 0.0;
-    elecstate::DensityMatrix<double, double> DM_trans(&this->kv, &this->pmat, this->nspin);
+    elecstate::DensityMatrix<double, double> DM_trans(&this->pmat, this->nspin, this->kv.kvec_d, this->nk);
     DM_trans.init_DMR(&GlobalC::GridD, &this->ucell);
     this->transition_dipole_.resize(X.get_nbands(), ModuleBase::Vector3<double>(0.0, 0.0, 0.0));
     for (int istate = 0;istate < X.get_nbands();++istate)
@@ -90,9 +90,9 @@ void LR::LR_Spectrum<std::complex<double>>::oscillator_strength()
     osc.resize(X.get_nbands(), 0.0);
     // const int nspin0 = (this->nspin == 2) ? 2 : 1;   use this in NSPIN=4 implementation
     double osc_tot = 0.0;
-    elecstate::DensityMatrix<std::complex<double>, std::complex<double>> DM_trans(&this->kv, &this->pmat, this->nspin);
+    elecstate::DensityMatrix<std::complex<double>, std::complex<double>> DM_trans(&this->pmat, this->nspin, this->kv.kvec_d, this->nk);
     DM_trans.init_DMR(&GlobalC::GridD, &this->ucell);
-    elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&this->kv, &this->pmat, this->nspin);
+    elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&this->pmat, this->nspin, this->kv.kvec_d, this->nk);
     DM_trans_real_imag.init_DMR(&GlobalC::GridD, &this->ucell);
 
     this->transition_dipole_.resize(X.get_nbands(), ModuleBase::Vector3<std::complex<double>>(0.0, 0.0, 0.0));
diff --git a/source/module_lr/lr_spectrum.h b/source/module_lr/lr_spectrum.h
index 4888e5e6b3..673a6b46ef 100644
--- a/source/module_lr/lr_spectrum.h
+++ b/source/module_lr/lr_spectrum.h
@@ -1,3 +1,4 @@
+#include "module_cell/klist.h"
 #include "module_lr/utils/gint_template.h"
 #include "module_psi/psi.h"
 #include "module_elecstate/module_dm/density_matrix.h"
diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp
index f9bc15f2f1..d3fffff390 100644
--- a/source/module_lr/operator_casida/operator_lr_hxc.cpp
+++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp
@@ -131,7 +131,7 @@ namespace LR
         ModuleBase::TITLE("OperatorLRHxc", "grid_calculation(complex)");
         ModuleBase::timer::tick("OperatorLRHxc", "grid_calculation");
 
-        elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&kv, pmat, nspin);
+        elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(pmat, nspin, kv.kvec_d, kv.get_nks() / nspin);
         DM_trans_real_imag.init_DMR(*this->hR);
         hamilt::HContainer<double> HR_real_imag(GlobalC::ucell, this->pmat);
         this->initialize_HR(HR_real_imag, ucell, gd, this->pmat);
diff --git a/source/module_lr/operator_casida/operator_lr_hxc.h b/source/module_lr/operator_casida/operator_lr_hxc.h
index 1c5ac17591..7c64f6c3df 100644
--- a/source/module_lr/operator_casida/operator_lr_hxc.h
+++ b/source/module_lr/operator_casida/operator_lr_hxc.h
@@ -1,4 +1,5 @@
 #pragma once
+#include "module_cell/klist.h"
 #include "module_hamilt_general/operator.h"
 #include "module_lr/utils/gint_template.h"
 #include "module_hamilt_lcao/module_gint/grid_technique.h"
@@ -87,7 +88,7 @@ namespace LR
                 for (int ib = prev_size;ib < nbands;++ib)
                 {
                     // the first dimenstion of DensityMatrix is nk=nks/nspin 
-                    DM_trans[ib] = LR_Util::make_unique<elecstate::DensityMatrix<T, TR>>(&this->kv, this->pmat, this->nspin);
+                    DM_trans[ib] = LR_Util::make_unique<elecstate::DensityMatrix<T, TR>>(this->pmat, this->nspin, this->kv.kvec_d, this->kv.get_nks() / nspin);
                     DM_trans[ib]->init_DMR(*this->hR);
                 }
             }

From 8f31794fb5bfe0e181094ba05f59b483e6d8f484 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Sun, 13 Oct 2024 14:52:20 +0800
Subject: [PATCH 2/6] rename _nks as _nk to avoid ambiguity

---
 .../module_dm/density_matrix.cpp              | 44 +++++++++----------
 .../module_dm/density_matrix.h                | 10 ++---
 source/module_io/td_current_io.cpp            |  2 +-
 source/module_lr/dm_trans/dmr_complex.cpp     |  4 +-
 4 files changed, 30 insertions(+), 30 deletions(-)

diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp
index 8ff877a741..3e2a1a9e93 100644
--- a/source/module_elecstate/module_dm/density_matrix.cpp
+++ b/source/module_elecstate/module_dm/density_matrix.cpp
@@ -26,10 +26,10 @@ DensityMatrix<TK, TR>::~DensityMatrix()
 
 template <typename TK, typename TR>
 DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin, const std::vector<ModuleBase::Vector3<double>>& kvec_d, const int nk)
-    : _paraV(paraV_in), _nspin(nspin), _kvec_d(kvec_d), _nks((nk > 0 && nk <= _kvec_d.size()) ? nk : _kvec_d.size())
+    : _paraV(paraV_in), _nspin(nspin), _kvec_d(kvec_d), _nk((nk > 0 && nk <= _kvec_d.size()) ? nk : _kvec_d.size())
 {
     ModuleBase::TITLE("DensityMatrix", "DensityMatrix-MK");
-    const int nks = _nks * _nspin;
+    const int nks = _nk * _nspin;
     this->_DMK.resize(nks);
     for (int ik = 0; ik < nks; ik++)
     {
@@ -39,7 +39,7 @@ DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const in
 }
 
 template <typename TK, typename TR>
-DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin) :_paraV(paraV_in), _nspin(nspin), _kvec_d({ ModuleBase::Vector3<double>(0,0,0) }), _nks(1)
+DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin) :_paraV(paraV_in), _nspin(nspin), _kvec_d({ ModuleBase::Vector3<double>(0,0,0) }), _nk(1)
 {
     ModuleBase::TITLE("DensityMatrix", "DensityMatrix-GO");
     this->_DMK.resize(_nspin);
@@ -236,7 +236,7 @@ template <typename TK, typename TR>
 TK* DensityMatrix<TK, TR>::get_DMK_pointer(const int ik) const
 {
 #ifdef __DEBUG
-    assert(ik < this->_nks * this->_nspin);
+    assert(ik < this->_nk * this->_nspin);
 #endif
     return const_cast<TK*>(this->_DMK[ik].data());
 }
@@ -246,7 +246,7 @@ template <typename TK, typename TR>
 void DensityMatrix<TK, TR>::set_DMK_pointer(const int ik, TK* DMK_in)
 {
 #ifdef __DEBUG
-    assert(ik < this->_nks * this->_nspin);
+    assert(ik < this->_nk * this->_nspin);
 #endif
     this->_DMK[ik].assign(DMK_in, DMK_in + this->_paraV->nrow * this->_paraV->ncol);
 }
@@ -257,17 +257,17 @@ void DensityMatrix<TK, TR>::set_DMK(const int ispin, const int ik, const int i,
 {
 #ifdef __DEBUG
     assert(ispin > 0 && ispin <= this->_nspin);
-    assert(ik >= 0 && ik < this->_nks);
+    assert(ik >= 0 && ik < this->_nk);
 #endif
     // consider transpose col=>row
-    this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->nrow + j] = value;
+    this->_DMK[ik + this->_nk * (ispin - 1)][i * this->_paraV->nrow + j] = value;
 }
 
 // set _DMK element
 template <typename TK, typename TR>
 void DensityMatrix<TK, TR>::set_DMK_zero()
 {
-    for (int ik = 0; ik < _nspin * _nks; ik++)
+    for (int ik = 0; ik < _nspin * _nk; ik++)
     {
         ModuleBase::GlobalFunc::ZEROS(this->_DMK[ik].data(),
                                       this->_paraV->get_row_size() * this->_paraV->get_col_size());
@@ -282,7 +282,7 @@ TK DensityMatrix<TK, TR>::get_DMK(const int ispin, const int ik, const int i, co
     assert(ispin > 0 && ispin <= this->_nspin);
 #endif
     // consider transpose col=>row
-    return this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->nrow + j];
+    return this->_DMK[ik + this->_nk * (ispin - 1)][i * this->_paraV->nrow + j];
 }
 
 // get _DMK nks, nrow, ncol
@@ -290,9 +290,9 @@ template <typename TK, typename TR>
 int DensityMatrix<TK, TR>::get_DMK_nks() const
 {
 #ifdef __DEBUG
-    assert(this->_DMK.size() == _nks * _nspin);
+    assert(this->_DMK.size() == _nk * _nspin);
 #endif
-    return _nks * _nspin;
+    return _nk * _nspin;
 }
 
 template <typename TK, typename TR>
@@ -364,7 +364,7 @@ void DensityMatrix<TK, TR>::cal_DMR_test()
 {
     for (int is = 1; is <= this->_nspin; ++is)
     {
-        int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2
+        int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
         hamilt::HContainer<TR>* tmp_DMR = this->_DMR[is - 1];
         // set zero since this function is called in every scf step
         tmp_DMR->set_zero();
@@ -396,7 +396,7 @@ void DensityMatrix<TK, TR>::cal_DMR_test()
 #endif
                 std::complex<TR> tmp_res;
                 // loop over k-points
-                for (int ik = 0; ik < this->_nks; ++ik)
+                for (int ik = 0; ik < this->_nk; ++ik)
                 {
                     // cal k_phase
                     // if TK==std::complex<double>, kphase is e^{ikR}
@@ -438,7 +438,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
     int ld_hk2 = 2 * ld_hk;
     for (int is = 1; is <= this->_nspin; ++is)
     {
-        int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2
+        int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
         hamilt::HContainer<double>* tmp_DMR = this->_DMR[is - 1];
         // set zero since this function is called in every scf step
         tmp_DMR->set_zero();
@@ -471,7 +471,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
                 // loop over k-points
                 if (PARAM.inp.nspin != 4)
                 {
-                    for (int ik = 0; ik < this->_nks; ++ik)
+                    for (int ik = 0; ik < this->_nk; ++ik)
                     {
                         // cal k_phase
                         // if TK==std::complex<double>, kphase is e^{ikR}
@@ -517,7 +517,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
                     std::vector<std::complex<double>> tmp_DMR(this->_paraV->get_col_size()
                                                                   * this->_paraV->get_row_size(),
                                                               std::complex<double>(0.0, 0.0));
-                    for (int ik = 0; ik < this->_nks; ++ik)
+                    for (int ik = 0; ik < this->_nk; ++ik)
                     {
                         // cal k_phase
                         // if TK==std::complex<double>, kphase is e^{ikR}
@@ -603,7 +603,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik)
     int ld_hk2 = 2 * ld_hk;
     for (int is = 1; is <= this->_nspin; ++is)
     {
-        int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2
+        int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
         hamilt::HContainer<double>* tmp_DMR = this->_DMR[is - 1];
         // set zero since this function is called in every scf step
         tmp_DMR->set_zero();
@@ -693,14 +693,14 @@ void DensityMatrix<double, double>::cal_DMR()
     int ld_hk = this->_paraV->nrow;
     for (int is = 1; is <= this->_nspin; ++is)
     {
-        int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2
+        int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
         hamilt::HContainer<double>* tmp_DMR = this->_DMR[is - 1];
         // set zero since this function is called in every scf step
         tmp_DMR->set_zero();
 
 #ifdef __DEBUG
         // assert(tmp_DMR->is_gamma_only() == true);
-        assert(this->_nks == 1);
+        assert(this->_nk == 1);
 #endif
 #ifdef _OPENMP
 #pragma omp parallel for
@@ -828,7 +828,7 @@ void DensityMatrix<TK, TR>::read_DMK(const std::string directory, const int ispi
     {
         for (int j = 0; j < this->_paraV->ncol; ++j)
         {
-            ifs >> this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j];
+            ifs >> this->_DMK[ik + this->_nk * (ispin - 1)][i * this->_paraV->ncol + j];
         }
     }
     ifs.close();
@@ -865,7 +865,7 @@ void DensityMatrix<double, double>::write_DMK(const std::string directory, const
             {
                 ofs << "\n";
             }
-            ofs << " " << this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j];
+            ofs << " " << this->_DMK[ik + this->_nk * (ispin - 1)][i * this->_paraV->ncol + j];
         }
     }
 
@@ -902,7 +902,7 @@ void DensityMatrix<std::complex<double>, double>::write_DMK(const std::string di
             {
                 ofs << "\n";
             }
-            ofs << " " << this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j].real();
+            ofs << " " << this->_DMK[ik + this->_nk * (ispin - 1)][i * this->_paraV->ncol + j].real();
         }
     }
 
diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h
index 02f389099a..56ee5737f9 100644
--- a/source/module_elecstate/module_dm/density_matrix.h
+++ b/source/module_elecstate/module_dm/density_matrix.h
@@ -243,8 +243,8 @@ class DensityMatrix
 
     /**
      * @brief density matrix in k space, which is a vector[ik]
-     * DMK should be a [_nspin][_nks][i][j] matrix,
-     * whose size is _nspin * _nks * _paraV->get_nrow() * _paraV->get_ncol()
+     * DMK should be a [_nspin][_nk][i][j] matrix,
+     * whose size is _nspin * _nk * _paraV->get_nrow() * _paraV->get_ncol()
      */
     // std::vector<ModuleBase::ComplexMatrix> _DMK;
     std::vector<std::vector<TK>> _DMK;
@@ -268,10 +268,10 @@ class DensityMatrix
 
     /**
      * @brief real number of k-points
-     * _nks is not equal to _kv->get_nks() when spin-polarization is considered
-     * _nks = kv->_nks / nspin
+     * _nk is not equal to _kv->get_nks() when spin-polarization is considered
+     * _nk = kv->get_nks() / nspin when nspin=2
      */
-    int _nks = 0;   // comment by lunasea: I want to rename it as nk
+    int _nk = 0;
 
 
 };
diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp
index 50fae430c2..31042e73d2 100644
--- a/source/module_io/td_current_io.cpp
+++ b/source/module_io/td_current_io.cpp
@@ -28,7 +28,7 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix<std::complex<double>, double>
     int ld_hk = DM_real.get_paraV_pointer()->nrow;
     int ld_hk2 = 2 * ld_hk;
     // tmp for is
-    int ik_begin = DM_real.get_DMK_nks() / nspin * (is - 1); // jump this->_nks for spin_down if nspin==2
+    int ik_begin = DM_real.get_DMK_nks() / nspin * (is - 1); // jump nk for spin_down if nspin==2
 
     hamilt::HContainer<double>* tmp_DMR_real = DM_real.get_DMR_vector()[is - 1];
     hamilt::HContainer<double>* tmp_DMR_imag = DM_imag.get_DMR_vector()[is - 1];
diff --git a/source/module_lr/dm_trans/dmr_complex.cpp b/source/module_lr/dm_trans/dmr_complex.cpp
index 265e962546..8cba19cb49 100644
--- a/source/module_lr/dm_trans/dmr_complex.cpp
+++ b/source/module_lr/dm_trans/dmr_complex.cpp
@@ -11,7 +11,7 @@ namespace elecstate
         ModuleBase::timer::tick("DensityMatrix", "cal_DMR");
         for (int is = 1; is <= this->_nspin; ++is)
         {
-            int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2
+            int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
             hamilt::HContainer<std::complex<double>>* tmp_DMR = this->_DMR[is - 1];
             // set zero since this function is called in every scf step
             tmp_DMR->set_zero();
@@ -43,7 +43,7 @@ namespace elecstate
 #endif
                     // loop over k-points
                     if (PARAM.inp.nspin != 4)
-                        for (int ik = 0; ik < this->_nks; ++ik)
+                        for (int ik = 0; ik < this->_nk; ++ik)
                         {
                             // cal k_phase
                             // if TK==std::complex<double>, kphase is e^{ikR}

From f15c6e1decc00f9dd6bdf646c3469aacb3a6b631 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci-lite[bot]"
 <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>
Date: Sun, 13 Oct 2024 07:47:04 +0000
Subject: [PATCH 3/6] [pre-commit.ci lite] apply automatic fixes

---
 source/module_lr/dm_trans/dmr_complex.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/source/module_lr/dm_trans/dmr_complex.cpp b/source/module_lr/dm_trans/dmr_complex.cpp
index 8cba19cb49..7a400f3bcc 100644
--- a/source/module_lr/dm_trans/dmr_complex.cpp
+++ b/source/module_lr/dm_trans/dmr_complex.cpp
@@ -42,7 +42,7 @@ namespace elecstate
                     }
 #endif
                     // loop over k-points
-                    if (PARAM.inp.nspin != 4)
+                    if (PARAM.inp.nspin != 4) {
                         for (int ik = 0; ik < this->_nk; ++ik)
                         {
                             // cal k_phase
@@ -70,9 +70,11 @@ namespace elecstate
                                 tmp_DMR_pointer += this->_paraV->get_col_size(iat2);
                             }
                         }
+}
                     // treat DMR as pauli matrix when NSPIN=4
-                    if (PARAM.inp.nspin == 4)
+                    if (PARAM.inp.nspin == 4) {
                         throw std::runtime_error("complex DM(R) with NSPIN=4 is not implemented yet");
+}
                 }
             }
         }

From 1c2e28c4cb04a3fb7c1479cc2f5402ae6fcdfa76 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Mon, 14 Oct 2024 14:50:35 +0800
Subject: [PATCH 4/6] remove map

---
 source/module_elecstate/elecstate_lcao.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp
index fdcbac0e43..6e30414531 100644
--- a/source/module_elecstate/elecstate_lcao.cpp
+++ b/source/module_elecstate/elecstate_lcao.cpp
@@ -131,7 +131,7 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
 template <typename TK>
 void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
 {
-    const int nspin_dm = std::map<int, int>({ {1, 1}, {2, 2}, {4, 1} })[nspin];
+    const int nspin_dm = nspin == 2 ? 2 : 1;
     this->DM = new DensityMatrix<TK, double>(paraV, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
 }
 

From 68158c94debccd967518b22f8f971a51750c065c Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Tue, 15 Oct 2024 18:37:29 +0800
Subject: [PATCH 5/6] make nk non-optional

---
 source/module_elecstate/module_dm/density_matrix.h            | 4 ++--
 source/module_elecstate/module_dm/test/test_dm_R_init.cpp     | 4 ++--
 .../module_elecstate/module_dm/test/test_dm_constructor.cpp   | 2 +-
 source/module_elecstate/module_dm/test/test_dm_io.cpp         | 4 ++--
 4 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h
index 56ee5737f9..f712a8e4c3 100644
--- a/source/module_elecstate/module_dm/density_matrix.h
+++ b/source/module_elecstate/module_dm/density_matrix.h
@@ -49,9 +49,9 @@ class DensityMatrix
      *  (usually {nspin_global -> nspin_dm} = {1->1, 2->2, 4->1}, but sometimes 2->1 like in LR-TDDFT)
      * @param kvec_d direct coordinates of kpoints
      * @param nk number of k-points, not always equal to K_Vectors::get_nks()/nspin_dm.
-     * if remains default or large than kvec_d.size(), it will be set to kvec_d.size()
+     *               it will be set to kvec_d.size() if the value is invalid
      */
-    DensityMatrix(const Parallel_Orbitals* _paraV, const int nspin, const std::vector<ModuleBase::Vector3<double>>& kvec_d, const int nk = -1);
+    DensityMatrix(const Parallel_Orbitals* _paraV, const int nspin, const std::vector<ModuleBase::Vector3<double>>& kvec_d, const int nk);
 
     /**
      * @brief Constructor of class DensityMatrix for gamma-only calculation, where kvector is not required
diff --git a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp
index e32308d456..b917e0ac6c 100644
--- a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp
+++ b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp
@@ -118,7 +118,7 @@ TEST_F(DMTest, DMInit1)
     // construct DM
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, nks);
     // initialize this->_DMR
     Grid_Driver gd(0,0);
     DM.init_DMR(&gd, &ucell);
@@ -145,7 +145,7 @@ TEST_F(DMTest, DMInit2)
     // construct DM
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, nks);
     // initialize Record_adj using Grid_Driver
     Grid_Driver gd(0,0);
     Record_adj ra;
diff --git a/source/module_elecstate/module_dm/test/test_dm_constructor.cpp b/source/module_elecstate/module_dm/test/test_dm_constructor.cpp
index c2c9832c44..8b4a9c7a6b 100644
--- a/source/module_elecstate/module_dm/test/test_dm_constructor.cpp
+++ b/source/module_elecstate/module_dm/test/test_dm_constructor.cpp
@@ -127,7 +127,7 @@ TEST_F(DMTest, DMConstructor_nspin1)
     std::cout << "dim0: " << paraV->dim0 << "    dim1:" << paraV->dim1 << std::endl;
     std::cout << "nrow: " << paraV->nrow << "    ncol:" << paraV->ncol << std::endl;
     int nspin = 1;
-    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, nks);
     // compare
     EXPECT_EQ(DM.get_DMK_nks(), kv->get_nks());
     EXPECT_EQ(DM.get_DMK_nrow(), paraV->nrow);
diff --git a/source/module_elecstate/module_dm/test/test_dm_io.cpp b/source/module_elecstate/module_dm/test/test_dm_io.cpp
index 8622bbc649..901ea7bef8 100644
--- a/source/module_elecstate/module_dm/test/test_dm_io.cpp
+++ b/source/module_elecstate/module_dm/test/test_dm_io.cpp
@@ -142,7 +142,7 @@ TEST_F(DMTest, DMConstructor1)
     int nspin = 1;
     // construct DM
     std::cout << paraV->nrow << paraV->ncol << std::endl;
-    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d);
+    elecstate::DensityMatrix<double, double> DM(paraV, nspin, kv->kvec_d, kv->get_nks());
     // read DMK
     std::string directory = "./support/";
     for (int is = 1; is <= nspin; ++is)
@@ -162,7 +162,7 @@ TEST_F(DMTest, DMConstructor1)
         }
     }
     // construct a new DM
-    elecstate::DensityMatrix<double, double> DM1(paraV, nspin, kv->kvec_d);
+    elecstate::DensityMatrix<double, double> DM1(paraV, nspin, kv->kvec_d, kv->get_nks());
     directory = "./support/output";
     for (int is = 1; is <= nspin; ++is)
     {

From 218e8ad31e426b21fd91c51e3bb3d5e40a380778 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Wed, 16 Oct 2024 19:02:44 +0800
Subject: [PATCH 6/6] fix after rebase

---
 source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp
index aacaa742dc..c8b9e01b94 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp
@@ -58,7 +58,8 @@ elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<d
 {
 
     // construct a DensityMatrix object
-    elecstate::DensityMatrix<std::complex<double>, double> edm(&kv, &pv, nspin);
+    const int nspin_dm = nspin == 2 ? 2 : 1;
+    elecstate::DensityMatrix<std::complex<double>, double> edm(&pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm);
 
     //--------------------------------------------
     // calculate the energy density matrix here.