fix: nspin = 2 for init_chg=dm

cmake/FindMKL.cmake

@@ -12,6 +12,10 @@ find_path(MKL_INCLUDE mkl_service.h HINTS ${MKLROOT}/include)
 find_library(MKL_INTEL NAMES mkl_intel_lp64 HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64)
 find_library(MKL_INTEL_THREAD NAMES mkl_intel_thread HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64)
 find_library(MKL_CORE NAMES mkl_core HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64)
+find_library(MKL_IOMP5 NAMES iomp5
+  HINTS ENV CMPLR_ROOT
+  PATH_SUFFIXES lib lib/intel64 linux/compiler/lib/intel64_lin
+)
 if(ENABLE_MPI)
   find_library(MKL_SCALAPACK NAMES mkl_scalapack_lp64 HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64)
   find_library(MKL_BLACS_INTELMPI NAMES mkl_blacs_intelmpi_lp64 HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64)
@@ -58,6 +62,11 @@ if(MKL_FOUND)
       IMPORTED_LOCATION "${MKL_BLACS_INTELMPI}"
       INTERFACE_INCLUDE_DIRECTORIES "${MKL_INCLUDE}")
   endif()
+  if(MKL_IOMP5 AND NOT TARGET MKL::IOMP5)
+    add_library(MKL::IOMP5 UNKNOWN IMPORTED)
+    set_target_properties(MKL::IOMP5 PROPERTIES
+      IMPORTED_LOCATION "${MKL_IOMP5}")
+  endif()
   add_library(MKL::MKL INTERFACE IMPORTED)
   if (ENABLE_MPI)
     set_property(TARGET MKL::MKL PROPERTY
@@ -74,6 +83,10 @@ if(MKL_FOUND)
     "-Wl,--end-group"
     )
   endif()
+  if(TARGET MKL::IOMP5)
+    set_property(TARGET MKL::MKL APPEND PROPERTY
+      INTERFACE_LINK_LIBRARIES MKL::IOMP5)
+  endif()
 endif()
 
 if(ENABLE_MPI)

docs/advanced/elec_properties/density_matrix.md

@@ -38,3 +38,15 @@ The following line is dimension of the density matrix, and the rest lines are th
 The examples can be found in [examples/density_matrix](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/density_matrix)
 
 - Note: now this function is valid only for LCAO gamma only calcualtion.
+
+## Real-space Density Matrix (CSR format)
+
+ABACUS can also output the real-space density matrix DM(R) in CSR (Compressed Sparse Row) format by setting:
+```
+out_dmr    1
+```
+After the calculation, the density matrix files are written to `OUT.${suffix}/`:
+- `nspin=1`: `dmrs1_nao.csr`
+- `nspin=2` (spin-polarized): `dmrs1_nao.csr` (spin-up) and `dmrs2_nao.csr` (spin-down)
+
+These files can be used to restart calculations by setting `init_chg dm` in the INPUT file together with `read_file_dir` pointing to the directory containing the CSR files. This is supported for both `nspin=1` and `nspin=2`.

docs/advanced/scf/initialization.md

@@ -1,15 +1,17 @@
 # Initializing SCF
 Good initializing would abate the number of iteration steps in SCF.
-Charge density should be initialed for constructing the initial hamiltonian operator. 
+Charge density should be initialed for constructing the initial hamiltonian operator.
 
 In PW basis, wavefunction should be initialized for iterate diagonalization method.
 In LCAO basis, wavefunction can be read to calculate initial charge density. The wavefunction itself does not have to be initialized.
 
 ## Charge Density
 `init_chg` is used for choosing the method of charge density initialization.
- - `atomic` : initial charge density by atomic charge density from pseudopotential file under keyword `PP_RHOATOM` 
+ - `atomic` : initial charge density by atomic charge density from pseudopotential file under keyword `PP_RHOATOM`
  - `file` : initial charge density from files produced by previous calculations with [`out_chg 1`](../elec_properties/charge.md).
  - `auto`: Abacus first attempts to read the density from a file; if not found, it defaults to using atomic density.
+ - `dm` (LCAO only): initial charge density from density matrix files in CSR format. For `nspin=1`, reads `dmrs1_nao.csr`. For `nspin=2` (spin-polarized), reads both `dmrs1_nao.csr` (spin-up) and `dmrs2_nao.csr` (spin-down). These files are generated by previous calculations with [`out_dmr 1`](../elec_properties/density_matrix.md). This method is particularly useful for restarting spin-polarized calculations.
+ - `hr` (LCAO only): initial charge density from Hamiltonian matrix files in CSR format. The Hamiltonian is read from file, then diagonalized to obtain wavefunctions and charge density. For `nspin=1`, reads `hrs1_nao.csr`. For `nspin=2` (spin-polarized), reads both `hrs1_nao.csr` (spin-up) and `hrs2_nao.csr` (spin-down). These files are generated by previous calculations with [`out_mat_hs2 1`](../input_files/input-main.md).
 
 ## Wave function
 `init_wfc` is used for choosing the method of wavefunction coefficient initialization.

source/source_esolver/esolver_ks_lcao.cpp

@@ -170,31 +170,21 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
 
     // 13.1) decide the strategy for initializing DMR and HR
-    if(istep == 0)//if the first scf step, readin DMR from file, 
+    if(istep == 0)//if the first scf step, readin DMR from file,
     {
         //calculate or readin the density matrix DMR
         if(PARAM.inp.init_chg == "dm")
         {
-            //! 13.1.1) init density matrix from file
-            std::string dmfile = PARAM.globalv.global_readin_dir + "dmrs1_nao.csr";
-            LCAO_domain::init_dm_from_file<TK>(dmfile, this->dmat, ucell, &(this->pv));
-	    GlobalV::ofs_running << " Read density matrix (real space) from file: " 
-		    << dmfile << std::endl;
+            //! 13.1.1) init charge density from density matrix file
+            LCAO_domain::init_chg_dm<TK>(PARAM.globalv.global_readin_dir, PARAM.inp.nspin,
+                this->dmat, ucell, &(this->pv), this->pelec->charge);
         }
         if(PARAM.inp.init_chg == "hr")
         {
-            //! 13.1.2) init HR from file
-            std::string hrfile = PARAM.globalv.global_readin_dir + "hrs1_nao.csr";
-            LCAO_domain::init_hr_from_file<TR>(
-                hrfile, 
-                dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt)->getHR(), 
-                ucell, 
-                &(this->pv)
-            );
-            this->p_hamilt->refresh(false);
-            hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
-            hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm, 
-                this->chr, PARAM.inp.nspin, 0);
+            //! 13.1.2) init charge density from Hamiltonian matrix file
+            LCAO_domain::init_chg_hr<TK, TR>(PARAM.globalv.global_readin_dir, PARAM.inp.nspin,
+                this->p_hamilt, ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm,
+                this->chr, PARAM.inp.ks_solver);
         }
     }
     else //if not, use the DMR calculated from last step
@@ -204,11 +194,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
         this->dmat.dm->cal_DMR();
     }
-    // 13.2 if init_chg = "dm", then calculate rho from readin DMR before init_scf
-    if(PARAM.inp.init_chg == "dm")
-    {
-        LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, this->pelec->charge, true);
-    }
     // 13.2) init_scf, should be before_scf? mohan add 2025-03-10
     this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 

source/source_lcao/CMakeLists.txt

@@ -64,6 +64,10 @@ if(ENABLE_LCAO)
     add_coverage(hamilt_lcao)
     endif()
 
+    IF (BUILD_TESTING)
+        add_subdirectory(test)
+    endif()
+
 endif()
 
 

source/source_lcao/LCAO_set.cpp

@@ -4,6 +4,9 @@
 #include "source_io/module_wf/read_wfc_nao.h" // use read_wfc_nao
 #include "source_estate/elecstate_tools.h" // use fixed_weights
 #include "source_lcao/module_hcontainer/read_hcontainer.h"
+#include "source_lcao/rho_tau_lcao.h" // use dm2rho
+#include "source_lcao/hamilt_lcao.h" // use HamiltLCAO for init_chg_hr
+#include "source_hsolver/hsolver_lcao.h" // use HSolverLCAO for init_chg_hr
 
 template <typename TK>
 void LCAO_domain::set_psi_occ_dm_chg(
@@ -94,20 +97,46 @@ void LCAO_domain::set_pot(
 
 template <typename TK>
 void LCAO_domain::init_dm_from_file(
-    const std::string dmfile,
+    const std::string& readin_dir,
+    const int nspin,
     LCAO_domain::Setup_DM<TK>& dmat,
     const UnitCell& ucell,
     const Parallel_Orbitals* pv)
 {
     ModuleBase::TITLE("LCAO_domain", "init_dm_from_file");
-    hamilt::HContainer<double>* dm_container = dmat.dm->get_DMR_vector()[0];
-    hamilt::Read_HContainer<double> reader_dm(
-        dm_container,
-        dmfile,
-        PARAM.globalv.nlocal,
-        &ucell
-    );
-    reader_dm.read();
+    const int nspin_dm = (nspin == 2) ? 2 : 1;
+    for (int is = 0; is < nspin_dm; ++is)
+    {
+        const std::string dmfile = readin_dir + "/dmrs" + std::to_string(is + 1) + "_nao.csr";
+        hamilt::HContainer<double>* dm_container = dmat.dm->get_DMR_vector()[is];
+        hamilt::Read_HContainer<double> reader_dm(
+            dm_container,
+            dmfile,
+            PARAM.globalv.nlocal,
+            &ucell
+        );
+        reader_dm.read();
+    }
+    return;
+}
+
+template <typename TK>
+void LCAO_domain::init_chg_dm(
+    const std::string& readin_dir,
+    const int nspin,
+    LCAO_domain::Setup_DM<TK>& dmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    Charge* chr)
+{
+    ModuleBase::TITLE("LCAO_domain", "init_chg_dm");
+
+    // Step 1: Read density matrix from file
+    LCAO_domain::init_dm_from_file<TK>(readin_dir, nspin, dmat, ucell, pv);
+
+    // Step 2: Convert density matrix to charge density
+    LCAO_domain::dm2rho(dmat.dm->get_DMR_vector(), nspin, chr, true);
+
     return;
 }
 
@@ -130,6 +159,57 @@ void LCAO_domain::init_hr_from_file(
     return;
 }
 
+template <typename TK, typename TR>
+void LCAO_domain::init_chg_hr(
+    const std::string& readin_dir,
+    const int nspin,
+    hamilt::Hamilt<TK>* p_hamilt,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    psi::Psi<TK>& psi,
+    elecstate::ElecState* pelec,
+    elecstate::DensityMatrix<TK, double>& dm,
+    Charge& chr,
+    const std::string& ks_solver)
+{
+    ModuleBase::TITLE("LCAO_domain", "init_chg_hr");
+
+    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(p_hamilt);
+    if (!hamilt_lcao)
+    {
+        ModuleBase::WARNING_QUIT("LCAO_domain::init_chg_hr", "p_hamilt is not HamiltLCAO");
+    }
+
+    // Step 1: Read HR from file(s)
+    if (nspin == 2)
+    {
+        // nspin=2: load spin-up into first half of hRS2, spin-down into second half
+        const std::string hrfile_up = readin_dir + "/hrs1_nao.csr";
+        LCAO_domain::init_hr_from_file<TR>(hrfile_up, hamilt_lcao->getHR(), ucell, pv);
+
+        // switch hR data pointer to spin-down half, then read hrs2
+        auto& hRS2 = hamilt_lcao->getHRS2();
+        hamilt_lcao->getHR()->allocate(hRS2.data() + hRS2.size() / 2, 0);
+        const std::string hrfile_down = readin_dir + "/hrs2_nao.csr";
+        LCAO_domain::init_hr_from_file<TR>(hrfile_down, hamilt_lcao->getHR(), ucell, pv);
+
+        // restore hR to spin-up half
+        hamilt_lcao->getHR()->allocate(hRS2.data(), 0);
+    }
+    else
+    {
+        const std::string hrfile = readin_dir + "/hrs1_nao.csr";
+        LCAO_domain::init_hr_from_file<TR>(hrfile, hamilt_lcao->getHR(), ucell, pv);
+    }
+
+    // Step 2: Mark HR as loaded from file (skip operator recalculation)
+    p_hamilt->refresh(false);
+
+    // Step 3: Diagonalize to get wavefunctions and charge density
+    hsolver::HSolverLCAO<TK> hsolver_lcao_obj(pv, ks_solver);
+    hsolver_lcao_obj.solve(p_hamilt, psi, pelec, dm, chr, nspin, 0);
+}
+
 
 
 template void LCAO_domain::set_psi_occ_dm_chg<double>(
@@ -183,16 +263,33 @@ template void LCAO_domain::set_pot<std::complex<double>>(
         const Input_para &inp);
 
 template void LCAO_domain::init_dm_from_file<double>(
-    const std::string dmfile,
+    const std::string& readin_dir,
+    const int nspin,
     LCAO_domain::Setup_DM<double>& dmat,
     const UnitCell& ucell,
     const Parallel_Orbitals* pv);
 template void LCAO_domain::init_dm_from_file<std::complex<double>>(
-    const std::string dmfile,
+    const std::string& readin_dir,
+    const int nspin,
     LCAO_domain::Setup_DM<std::complex<double>>& dmat,
     const UnitCell& ucell,
     const Parallel_Orbitals* pv);
 
+template void LCAO_domain::init_chg_dm<double>(
+    const std::string& readin_dir,
+    const int nspin,
+    LCAO_domain::Setup_DM<double>& dmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    Charge* chr);
+template void LCAO_domain::init_chg_dm<std::complex<double>>(
+    const std::string& readin_dir,
+    const int nspin,
+    LCAO_domain::Setup_DM<std::complex<double>>& dmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    Charge* chr);
+
 template void LCAO_domain::init_hr_from_file<double>(
     const std::string hrfile,
     hamilt::HContainer<double>* hmat,
@@ -203,3 +300,37 @@ template void LCAO_domain::init_hr_from_file<std::complex<double>>(
     hamilt::HContainer<std::complex<double>>* hmat,
     const UnitCell& ucell,
     const Parallel_Orbitals* pv);
+
+template void LCAO_domain::init_chg_hr<double, double>(
+    const std::string& readin_dir,
+    const int nspin,
+    hamilt::Hamilt<double>* p_hamilt,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    psi::Psi<double>& psi,
+    elecstate::ElecState* pelec,
+    elecstate::DensityMatrix<double, double>& dm,
+    Charge& chr,
+    const std::string& ks_solver);
+template void LCAO_domain::init_chg_hr<std::complex<double>, double>(
+    const std::string& readin_dir,
+    const int nspin,
+    hamilt::Hamilt<std::complex<double>>* p_hamilt,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    psi::Psi<std::complex<double>>& psi,
+    elecstate::ElecState* pelec,
+    elecstate::DensityMatrix<std::complex<double>, double>& dm,
+    Charge& chr,
+    const std::string& ks_solver);
+template void LCAO_domain::init_chg_hr<std::complex<double>, std::complex<double>>(
+    const std::string& readin_dir,
+    const int nspin,
+    hamilt::Hamilt<std::complex<double>>* p_hamilt,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv,
+    psi::Psi<std::complex<double>>& psi,
+    elecstate::ElecState* pelec,
+    elecstate::DensityMatrix<std::complex<double>, double>& dm,
+    Charge& chr,
+    const std::string& ks_solver);