Feature: support init_chg dm to restart with DMR

source/Makefile.Objects

@@ -356,6 +356,7 @@ OBJS_HCONTAINER=base_matrix.o\
     atom_pair.o\
     hcontainer.o\
     output_hcontainer.o\
+    read_hcontainer.o\
     func_folding.o\
     func_transfer.o\
     transfer.o\

source/source_esolver/esolver_ks_lcao.cpp

@@ -75,6 +75,13 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
 
     LCAO_domain::set_psi_occ_dm_chg<TK>(this->kv, this->psi, this->pv, this->pelec,
       this->dmat, this->chr, inp);
+
+    if(inp.init_chg == "dm")
+    {
+        //! 4.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));
+    }
 
     LCAO_domain::set_pot<TK>(ucell, this->kv, this->sf, *this->pw_rho, *this->pw_rhod,
       this->pelec, this->orb_, this->pv, this->locpp, this->dftu,
@@ -160,7 +167,12 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
     this->exx_nao.before_scf(ucell, this->kv, orb_, this->p_chgmix, istep, PARAM.inp);
 
-    // 12) init_scf, should be before_scf? mohan add 2025-03-10
+    // 12.1) 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);
+    }
+    // 12.2) init_scf, should be before_scf? mohan add 2025-03-10
     this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
     // 13) initalize DM(R), which has the same size with Hamiltonian(R)
@@ -169,7 +181,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     {
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","p_hamilt does not exist");
     }
-    this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
+    if(PARAM.inp.init_chg != "dm") this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
 
 #ifdef __MLALGO
     // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R)

source/source_io/read_input_item_system.cpp

@@ -552,7 +552,7 @@ void ReadInput::item_system()
             }
         };
         item.check_value = [](const Input_Item& item, const Parameter& para) {
-            const std::vector<std::string> init_chgs = {"atomic", "file", "wfc", "auto"};
+            const std::vector<std::string> init_chgs = {"atomic", "file", "wfc", "auto", "dm"};
             if (std::find(init_chgs.begin(), init_chgs.end(), para.input.init_chg) == init_chgs.end())
             {
                 const std::string warningstr = nofound_str(init_chgs, "init_chg");

source/source_lcao/LCAO_set.cpp

@@ -3,6 +3,7 @@
 #include "source_psi/setup_psi.h" // use Setup_Psi
 #include "source_io/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"
 
 template <typename TK>
 void LCAO_domain::set_psi_occ_dm_chg(
@@ -91,6 +92,27 @@ void LCAO_domain::set_pot(
     return;
 }
 
+template <typename TK>
+void LCAO_domain::init_dm_from_file(
+    const std::string dmfile,
+    LCAO_domain::Setup_DM<TK>& dmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv)
+{
+    ModuleBase::TITLE("LCAO_domain::init_dm_from_file", "init_dm_from_file");
+    hamilt::HContainer<double>* dm_container = new hamilt::HContainer<double>(pv);
+    dmat.dm->init_DMR(dm_container[0]);
+    hamilt::Read_HContainer<double> reader_dm(
+        dmat.dm->get_DMR_vector()[0],
+        dmfile,
+        PARAM.globalv.nlocal,
+        &ucell
+    );
+    reader_dm.read();
+    delete dm_container;
+    return;
+}
+
 
 
 template void LCAO_domain::set_psi_occ_dm_chg<double>(
@@ -142,3 +164,14 @@ template void LCAO_domain::set_pot<std::complex<double>>(
         Exx_NAO<std::complex<double>> &exx_nao,
         Setup_DeePKS<std::complex<double>> &deepks,
         const Input_para &inp);
+
+template void LCAO_domain::init_dm_from_file<double>(
+    const std::string dmfile,
+    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,
+    LCAO_domain::Setup_DM<std::complex<double>>& dmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv);

source/source_lcao/LCAO_set.h

@@ -53,6 +53,12 @@ void set_pot(
         Setup_DeePKS<TK> &deepks,
         const Input_para &inp);
 
+template <typename TK>
+void init_dm_from_file(
+	const std::string dmfile,
+	LCAO_domain::Setup_DM<TK>& dmat,
+	const UnitCell& ucell,
+	const Parallel_Orbitals* pv);
 } // end namespace
 
 #endif

source/source_lcao/module_hcontainer/CMakeLists.txt

@@ -5,6 +5,7 @@ list(APPEND objects
     atom_pair.cpp
     hcontainer.cpp
     output_hcontainer.cpp
+    read_hcontainer.cpp
     func_folding.cpp
     transfer.cpp
     func_transfer.cpp

source/source_lcao/module_hcontainer/read_hcontainer.cpp

@@ -0,0 +1,169 @@
+#include "read_hcontainer.h"
+
+#include "source_io/sparse_matrix.h"
+#include "source_io/csr_reader.h"
+#include "hcontainer_funcs.h"
+
+#include <fstream>
+
+namespace hamilt
+{
+
+/**
+ * @brief Constructor of Read_HContainer
+ * @attention ifs should be open outside of this interface
+ */
+template <typename T>
+Read_HContainer<T>::Read_HContainer(hamilt::HContainer<T>* hcontainer,
+                                        const std::string& filename,
+                                        const int nlocal,
+                                        const UnitCell* ucell)
+    : _hcontainer(hcontainer), _filename(filename), _nlocal(nlocal), _ucell(ucell)
+{
+}
+
+template <typename T>
+void Read_HContainer<T>::read()
+{
+    // build atom index of col and row
+    std::vector<int> atom_index_row;
+    std::vector<int> atom_index_col;
+    int natom = this->_ucell->nat;
+    Parallel_Orbitals pv_serial;
+    pv_serial.set_serial(this->_nlocal, this->_nlocal);
+    pv_serial.set_atomic_trace(this->_ucell->get_iat2iwt(), this->_ucell->nat, this->_nlocal);
+    for (int iat = 0; iat < natom; ++iat)
+    {
+        int row_size = pv_serial.get_row_size(iat);
+        int col_size = pv_serial.get_col_size(iat);
+        for (int i = 0; i < row_size; ++i)
+        {
+            atom_index_row.push_back(iat);
+        }
+        for (int j = 0; j < col_size; ++j)
+        {
+            atom_index_col.push_back(iat);
+        }
+    }
+    //
+    hamilt::HContainer<T> hcontainer_serial(&pv_serial);
+
+#ifdef __MPI
+    if(GlobalV::MY_RANK == 0)
+    {
+#endif
+    ModuleIO::csrFileReader<T> csr(this->_filename);
+    int step = csr.getStep();
+    int matrix_dimension = csr.getMatrixDimension();
+    int r_number = csr.getNumberOfR();
+
+    //construct serial hcontainer firstly
+    // prepare atom index mapping from csr row/col to atom index
+    for (int i = 0; i < r_number; i++)
+    {
+        std::vector<int> RCoord = csr.getRCoordinate(i);
+        ModuleIO::SparseMatrix<T> sparse_matrix = csr.getMatrix(i);
+        for (const auto& element: sparse_matrix.getElements())
+        {
+            int row = element.first.first;
+            int col = element.first.second;
+            T value = element.second;
+
+
+            //insert into hcontainer
+            int atom_i = atom_index_row[row];
+            int atom_j = atom_index_col[col];
+            auto* ij_pair = hcontainer_serial.find_pair(atom_i, atom_j);
+            if(ij_pair == nullptr)
+            {
+                //insert new pair
+                hamilt::AtomPair<T> new_pair(atom_i, atom_j, RCoord[0], RCoord[1], RCoord[2], &pv_serial);
+                hcontainer_serial.insert_pair(new_pair);
+            }
+            else
+            {
+                if(ij_pair->find_R(RCoord[0], RCoord[1], RCoord[2]) == -1)
+                {
+                    //insert new R
+                    hamilt::AtomPair<T> new_pair(atom_i, atom_j, RCoord[0], RCoord[1], RCoord[2], &pv_serial);
+                    hcontainer_serial.insert_pair(new_pair);
+                }
+            }
+        }
+    }
+    hcontainer_serial.allocate(nullptr, true);
+    // second loop, add values into hcontainer
+    for (int i = 0; i < r_number; i++)
+    {
+        std::vector<int> RCoord = csr.getRCoordinate(i);
+        ModuleIO::SparseMatrix<T> sparse_matrix = csr.getMatrix(i);
+        for (const auto& element: sparse_matrix.getElements())
+        {
+            int row = element.first.first;
+            int col = element.first.second;
+            T value = element.second;
+
+            //insert into hcontainer
+            int atom_i = atom_index_row[row];
+            int atom_j = atom_index_col[col];
+            auto* matrix = hcontainer_serial.find_matrix(atom_i, atom_j, RCoord[0], RCoord[1], RCoord[2]);
+            matrix->add_element(row - pv_serial.atom_begin_row[atom_i],
+                                col - pv_serial.atom_begin_col[atom_j],
+                                value);
+        }
+    }
+#ifdef __MPI
+}
+    // thirdly, distribute hcontainer_serial to parallel hcontainer
+    // send <IJR>s from serial_rank to all ranks
+    int my_rank, size;
+    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    std::vector<int> para_ijrs;
+    if (my_rank == 0)
+    {
+        para_ijrs = hcontainer_serial.get_ijr_info();
+        this->_hcontainer->insert_ijrs(&para_ijrs);
+        this->_hcontainer->allocate();
+    }
+    if (my_rank != 0)
+    {
+        std::vector<int> tmp_ijrs;
+        MPI_Status status;
+        long tmp_size = 0;
+        MPI_Recv(&tmp_size, 1, MPI_LONG, 0, 0, MPI_COMM_WORLD, &status);
+        tmp_ijrs.resize(tmp_size);
+        MPI_Recv(tmp_ijrs.data(),
+                    tmp_ijrs.size(),
+                    MPI_INT,
+                    0,
+                    1,
+                    MPI_COMM_WORLD,
+                    &status);
+        this->_hcontainer->insert_ijrs(&tmp_ijrs);
+        this->_hcontainer->allocate();
+    }
+    else
+    {
+        for (int i = 1; i < size; ++i)
+        {
+            long tmp_size = para_ijrs.size();
+            MPI_Send(&tmp_size, 1, MPI_LONG, i, 0, MPI_COMM_WORLD);
+            MPI_Send(para_ijrs.data(), para_ijrs.size(), MPI_INT, i, 1, MPI_COMM_WORLD);
+        }
+    }
+    // gather values from serial_rank to Parallels
+    transferSerial2Parallels(hcontainer_serial, this->_hcontainer, 0);
+#else
+    std::vector<int> para_ijrs = hcontainer_serial.get_ijr_info();
+    this->_hcontainer->insert_ijrs(&para_ijrs);
+    this->_hcontainer->allocate();
+    this->_hcontainer->add(hcontainer_serial);
+#endif
+
+}
+
+template class Read_HContainer<double>;
+template class Read_HContainer<std::complex<double>>;
+
+} // namespace hamilt

source/source_lcao/module_hcontainer/read_hcontainer.h

@@ -0,0 +1,47 @@
+#ifndef READ_HCONTAINER_H
+#define READ_HCONTAINER_H
+
+#include "source_lcao/module_hcontainer/hcontainer.h"
+#include "source_cell/unitcell.h"
+
+namespace hamilt
+{
+
+/**
+ * @brief A class to read the HContainer
+ */
+template <typename T>
+class Read_HContainer
+{
+  public:
+    Read_HContainer(
+        hamilt::HContainer<T>* hcontainer, 
+        const std::string& filename,
+        const int nlocal,
+        const UnitCell* ucell
+    );
+    // read the matrices of all R vectors to the read stream
+    void read();
+
+    /**
+     * read the matrix of a single R vector to the output stream
+     * rx_in, ry_in, rz_in: the R vector from the input
+     */
+    void read(int rx_in, int ry_in, int rz_in);
+
+    /**
+     * read the matrix of a single R vector to the output stream
+     * rx, ry, rz: the R vector from the HContainer
+     */
+    void read_single_R(int rx, int ry, int rz);
+
+  private:
+    hamilt::HContainer<T>* _hcontainer;
+    std::string _filename;
+    int _nlocal;
+    const UnitCell* _ucell;
+};
+
+} // namespace hamilt
+
+#endif // OUTPUT_HCONTAINER_H

source/source_lcao/rho_tau_lcao.cpp

@@ -4,7 +4,8 @@
 
 void LCAO_domain::dm2rho(std::vector<hamilt::HContainer<double>*> &dmr,
     const int nspin,
-    Charge* chr)
+    Charge* chr,
+    bool skip_normalize)
 {
     ModuleBase::TITLE("LCAO_domain", "dm2rho");
     ModuleBase::timer::tick("LCAO_domain", "dm2rho");
@@ -16,7 +17,7 @@ void LCAO_domain::dm2rho(std::vector<hamilt::HContainer<double>*> &dmr,
 
     ModuleGint::cal_gint_rho(dmr, nspin, chr->rho);
 
-    chr->renormalize_rho();
+    if(!skip_normalize)chr->renormalize_rho();
 
     // should be moved somewhere else, mohan 20251024
 	if (XC_Functional::get_ked_flag())

source/source_lcao/rho_tau_lcao.h

@@ -9,7 +9,8 @@ namespace LCAO_domain
 {
 	void dm2rho(std::vector<hamilt::HContainer<double>*> &dmr,
 			const int nspin,
-			Charge* chr);
+			Charge* chr, 
+			bool skip_normalize = false);
 
 	void dm2tau(std::vector<hamilt::HContainer<double>*> &dmr,
 			const int nspin,

tests/03_NAO_multik/02_NO_KP_15/INPUT

@@ -4,7 +4,7 @@ suffix			autotest
 calculation     scf
 
 nbands			6
-symmetry		1
+symmetry		0
 pseudo_dir	../../PP_ORB
 orbital_dir	../../PP_ORB
 gamma_only      0
@@ -21,6 +21,9 @@ basis_type		lcao
 #Parameters (4.Smearing)
 smearing_method		gauss
 smearing_sigma			0.002
+init_chg                dm
+#out_dmr                  1
+read_file_dir           ./
 
 #Parameters (5.Mixing)
 mixing_type		broyden