Feature: init_chg hr for initializing charge density with hrs1_nao.csr (deepmodeling#6791)

* Fix: init_chg dm with mismatched sparsity

* Feature:  for initializing charge density with hrs1_nao.csr

* Test: add a test case

---------

Co-authored-by: dyzheng <zhengdy@bjaisi.com>

Author: dyzheng

docs/advanced/input_files/input-main.md

@@ -640,6 +640,8 @@ These variables are used to control general system parameters.
   - atomic: the density is starting from the summation of the atomic density of single atoms.
   - file: the density will be read in from a binary file `charge-density.dat` first. If it does not exist, the charge density will be read in from cube files. Besides, when you do `nspin=1` calculation, you only need the density file chgs1.cube. However, if you do `nspin=2` calculation, you also need the density file chgs2.cube. The density file should be output with these names if you set out_chg = 1 in INPUT file.
   - wfc: the density will be calculated by wavefunctions and occupations. Wavefunctions are read in from binary files `wf*.dat` (see [out_wfc_pw](#out_wfc_pw)) while occupations are read in from file `eig.txt`.
+  - dm: the density will be calculated by real space density matrix(DMR) of LCAO base. DMR is read in from file `dmrs1_nao.csr` in directory [read_file_dir](#read_file_dir).
+  - hr: the real space Hamiltonian matrix(HR) will be read in from file `hrs1_nao.csr` in directory [read_file_dir](#read_file_dir), and DMR and charge density will be calculated from it.
   - auto: Abacus first attempts to read the density from a file; if not found, it defaults to using atomic density.
 - **Default**: atomic
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -75,13 +75,6 @@ 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,
@@ -167,46 +160,69 @@ 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.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(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
-
-    // 13) initalize DM(R), which has the same size with Hamiltonian(R)
+    // 12) initalize DM(R), which has the same size with Hamiltonian(R)
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
+
     if(!hamilt_lcao)
     {
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","p_hamilt does not exist");
     }
-    if(PARAM.inp.init_chg != "dm") this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
+    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, 
+    {
+        //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));
+        }
+        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);
+        }
+    }
+    else //if not, use the DMR calculated from last step
+    {
+        // 13.1.2) two cases are considered:
+        // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
+        // 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);
 
 #ifdef __MLALGO
     // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R)
     this->deepks.ld.init_DMR(ucell, orb_, this->pv, this->gd);
 #endif
 
-    // 15) two cases are considered:
-    // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
-    // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
-    if (istep > 0)
-    {
-        this->dmat.dm->cal_DMR();
-    }
-
     // 16) the electron charge density should be symmetrized,
     Symmetry_rho srho;
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
         srho.begin(is, this->chr, this->pw_rho, ucell.symm);
     }
 
-    // 17) why we need to set this sentence? mohan add 2025-03-10
-    this->p_hamilt->non_first_scf = istep;
-
-    // 18) update of RDMFT, added by jghan
+    // 17) update of RDMFT, added by jghan
     if (PARAM.inp.rdmft == true)
     {
         rdmft_solver.update_ion(ucell, *(this->pw_rho), this->locpp.vloc, this->sf.strucFac);

source/source_hamilt/hamilt.h

@@ -21,7 +21,7 @@ class Hamilt
     virtual void updateHk(const int ik){return;}
 
     /// refresh status of Hamiltonian, for example, refresh H(R) and S(R) in LCAO case
-    virtual void refresh(void){return;}
+    virtual void refresh(bool yes = true){return;}
 
     /// core function: for solving eigenvalues of Hamiltonian with iterative method
 	virtual void hPsi(
@@ -55,8 +55,6 @@ class Hamilt
 
     std::string classname = "none";
 
-    int non_first_scf=0;
-
     /// first node operator, add operations from each operators
     Operator<T, Device>* ops = nullptr;
 

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", "dm"};
+            const std::vector<std::string> init_chgs = {"atomic", "file", "wfc", "auto", "dm", "hr"};
             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

@@ -99,17 +99,34 @@ void LCAO_domain::init_dm_from_file(
     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]);
+    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(
-        dmat.dm->get_DMR_vector()[0],
+        dm_container,
         dmfile,
         PARAM.globalv.nlocal,
         &ucell
     );
     reader_dm.read();
-    delete dm_container;
+    return;
+}
+
+template <typename TR>
+void LCAO_domain::init_hr_from_file(
+    const std::string hrfile,
+    hamilt::HContainer<TR>* hmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv)
+{
+    ModuleBase::TITLE("LCAO_domain", "init_hr_from_file");
+    hmat->set_zero();
+    hamilt::Read_HContainer<TR> reader_hr(
+        hmat,
+        hrfile,
+        PARAM.globalv.nlocal,
+        &ucell
+    );
+    reader_hr.read();
     return;
 }
 
@@ -175,3 +192,14 @@ template void LCAO_domain::init_dm_from_file<std::complex<double>>(
     LCAO_domain::Setup_DM<std::complex<double>>& dmat,
     const UnitCell& ucell,
     const Parallel_Orbitals* pv);
+
+template void LCAO_domain::init_hr_from_file<double>(
+    const std::string hrfile,
+    hamilt::HContainer<double>* hmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv);
+template void LCAO_domain::init_hr_from_file<std::complex<double>>(
+    const std::string hrfile,
+    hamilt::HContainer<std::complex<double>>* hmat,
+    const UnitCell& ucell,
+    const Parallel_Orbitals* pv);

source/source_lcao/LCAO_set.h

@@ -53,12 +53,25 @@ void set_pot(
         Setup_DeePKS<TK> &deepks,
         const Input_para &inp);
 
+/**
+ * @brief read in DMR from file, and save it into dmat
+ */
 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);
+
+/**
+ * @brief read in HR from file, and save it into hmat
+ */
+template <typename TK>
+void init_hr_from_file(
+	const std::string hrfile,
+	hamilt::HContainer<TK>* hmat,
+	const UnitCell& ucell,
+	const Parallel_Orbitals* pv);
 } // end namespace
 
 #endif

source/source_lcao/hamilt_lcao.cpp

@@ -493,20 +493,32 @@ void HamiltLCAO<TK, TR>::updateHk(const int ik)
 }
 
 template <typename TK, typename TR>
-void HamiltLCAO<TK, TR>::refresh()
+void HamiltLCAO<TK, TR>::refresh(bool yes)
 {
     ModuleBase::TITLE("HamiltLCAO", "refresh");
-    dynamic_cast<hamilt::OperatorLCAO<TK, TR>*>(this->ops)->set_hr_done(false);
-    if (PARAM.inp.nspin == 2)
+    if(yes)
     {
-        this->refresh_times = 1;
-        this->current_spin = 0;
-        if (this->hR->get_nnr() != this->hRS2.size() / 2)
+        dynamic_cast<hamilt::OperatorLCAO<TK, TR>*>(this->ops)->set_hr_done(false);
+        if (PARAM.inp.nspin == 2)
         {
-            // operator has changed, resize hRS2
-            this->hRS2.resize(this->hR->get_nnr() * 2);
+            this->refresh_times = 1;
+            this->current_spin = 0;
+            if (this->hR->get_nnr() != this->hRS2.size() / 2)
+            {
+                // operator has changed, resize hRS2
+                this->hRS2.resize(this->hR->get_nnr() * 2);
+            }
+            this->hR->allocate(this->hRS2.data(), 0);
+        }
+    }
+    else {
+        dynamic_cast<hamilt::OperatorLCAO<TK, TR>*>(this->ops)->set_hr_done(true);
+        this->refresh_times = 0;
+        if (PARAM.inp.nspin == 2)
+        {
+            ModuleBase::WARNING_QUIT("HamiltLCAO::refresh",
+                                      "When turning off the refresh flag, the nspin==2 case is not supported yet.");
         }
-        this->hR->allocate(this->hRS2.data(), 0);
     }
 }
 

source/source_lcao/hamilt_lcao.h

@@ -124,7 +124,7 @@ class HamiltLCAO : public Hamilt<TK>
 #endif
 
     /// refresh the status of HR
-    void refresh() override;
+    void refresh(bool yes) override;
 
     // for target K point, update consequence of hPsi() and matrix()
     virtual void updateHk(const int ik) override;

source/source_lcao/module_operator_lcao/operator_lcao.cpp

@@ -98,7 +98,7 @@ void OperatorLCAO<TK, TR>::init(const int ik_in) {
         // cal_type=lcao_overlap refer to overlap matrix operators, which are
         // only rely on stucture, and not changed during SCF
 
-        if (!this->hr_done) {
+        {
             // update SR first
             // in cal_type=lcao_overlap, SR should be updated by each sub-chain
             // nodes
@@ -228,7 +228,7 @@ void OperatorLCAO<TK, TR>::init(const int ik_in) {
         != nullptr) { // it is not the last node, loop next init() function
         // pass HR status to next node and than set HR status of this node to
         // done
-        if (!this->hr_done) {
+        {
             dynamic_cast<OperatorLCAO<TK, TR>*>(this->next_op)->hr_done
                 = this->hr_done;
         }

tests/03_NAO_multik/59_NO_KP_initHR/INPUT

@@ -0,0 +1,30 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix			autotest
+calculation     scf
+
+nbands			6
+symmetry		0
+
+#Parameters (2.Iteration)
+ecutwfc			20
+scf_thr				1e-6
+scf_nmax			100
+
+
+#Parameters (3.Basis)
+basis_type		lcao
+
+#Parameters (4.Smearing)
+smearing_method		gauss
+smearing_sigma			0.002
+
+#Parameters (5.Mixing)
+mixing_type		broyden
+mixing_beta		0.7
+pseudo_dir	../../PP_ORB
+orbital_dir	../../PP_ORB
+
+#out_mat_hs2     1
+init_chg        hr
+read_file_dir   ./