Refactor: set dm_to_rho as a new esolver

source/Makefile.Objects

@@ -266,6 +266,7 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       lcao_after_scf.o\
       esolver_gets.o\
       lcao_others.o\
+      esolver_dm2rho.o\
 
 OBJS_GINT=gint.o\
       gint_gamma_env.o\

source/module_elecstate/elecstate_lcao.cpp

@@ -23,32 +23,6 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d
     ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
     ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
 
-    // // the calculations of dm, and dm -> rho are, technically, two separate
-    // // functionalities, as we cannot rule out the possibility that we may have a
-    // // dm from other sources, such as read from file. However, since we are not
-    // // separating them now, I opt to add a flag to control how dm is obtained as
-    // // of now
-    // if (!PARAM.inp.dm_to_rho)
-    // {
-    //     ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
-
-    //     // this part for calculating DMK in 2d-block format, not used for charge
-    //     // now
-    //     //    psi::Psi<std::complex<double>> dm_k_2d();
-
-    //     if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver ==
-    //     "scalapack_gvx" || PARAM.inp.ks_solver == "lapack"
-    //         || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp"
-    //         || PARAM.inp.ks_solver == "cg_in_lcao") // Peize Lin test 2019-05-15
-    //     {
-    //         elecstate::cal_dm_psi(this->DM->get_paraV_pointer(),
-    //                               this->wg,
-    //                               psi,
-    //                               *(this->DM));
-    //         this->DM->cal_DMR();
-    //     }
-    // }
-
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
         ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],

source/module_esolver/CMakeLists.txt

@@ -20,6 +20,7 @@ if(ENABLE_LCAO)
       lcao_after_scf.cpp
       esolver_gets.cpp
       lcao_others.cpp
+      esolver_dm2rho.cpp
   )
 endif()
 

source/module_esolver/esolver.cpp

@@ -5,6 +5,7 @@
 #include "module_base/module_device/device.h"
 #include "module_parameter/parameter.h"
 #ifdef __LCAO
+#include "esolver_dm2rho.h"
 #include "esolver_gets.h"
 #include "esolver_ks_lcao.h"
 #include "esolver_ks_lcao_tddft.h"
@@ -204,11 +205,25 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
         }
         else if (PARAM.inp.nspin < 4)
         {
-            return new ESolver_KS_LCAO<std::complex<double>, double>();
+            if (PARAM.inp.dm_to_rho)
+            {
+                return new ESolver_DM2rho<std::complex<double>, double>();
+            }
+            else
+            {
+                return new ESolver_KS_LCAO<std::complex<double>, double>();
+            }
         }
         else
         {
-            return new ESolver_KS_LCAO<std::complex<double>, std::complex<double>>();
+            if (PARAM.inp.dm_to_rho)
+            {
+                return new ESolver_DM2rho<std::complex<double>, std::complex<double>>();
+            }
+            else
+            {
+                return new ESolver_KS_LCAO<std::complex<double>, std::complex<double>>();
+            }
         }
     }
     else if (esolver_type == "ksdft_lcao_tddft")

source/module_esolver/esolver_dm2rho.cpp

@@ -0,0 +1,100 @@
+#include "esolver_dm2rho.h"
+
+#include "module_base/timer.h"
+#include "module_cell/module_neighbor/sltk_atom_arrange.h"
+#include "module_elecstate/elecstate_lcao.h"
+#include "module_elecstate/read_pseudo.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
+#include "module_io/cube_io.h"
+#include "module_io/io_npz.h"
+#include "module_io/print_info.h"
+
+namespace ModuleESolver
+{
+
+template <typename TK, typename TR>
+ESolver_DM2rho<TK, TR>::ESolver_DM2rho()
+{
+    this->classname = "ESolver_DM2rho";
+    this->basisname = "LCAO";
+}
+
+template <typename TK, typename TR>
+ESolver_DM2rho<TK, TR>::~ESolver_DM2rho()
+{
+}
+
+template <typename TK, typename TR>
+void ESolver_DM2rho<TK, TR>::before_all_runners(UnitCell& ucell, const Input_para& inp)
+{
+    ModuleBase::TITLE("ESolver_DM2rho", "before_all_runners");
+    ModuleBase::timer::tick("ESolver_DM2rho", "before_all_runners");
+
+    ESolver_KS_LCAO<TK, TR>::before_all_runners(ucell, inp);
+
+    ModuleBase::timer::tick("ESolver_DM2rho", "before_all_runners");
+}
+
+template <typename TK, typename TR>
+void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
+{
+    ModuleBase::TITLE("ESolver_DM2rho", "runner");
+    ModuleBase::timer::tick("ESolver_DM2rho", "runner");
+
+    ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
+
+    // file name of DM
+    std::string zipname = "output_DM0.npz";
+    elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
+
+    // read DM from file
+    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
+
+    // if nspin=2, need extra reading
+    if (PARAM.inp.nspin == 2)
+    {
+        zipname = "output_DM1.npz";
+        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
+    }
+
+    this->pelec->psiToRho(*this->psi);
+
+    int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
+
+    for (int is = 0; is < nspin0; is++)
+    {
+        std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
+
+        // write electron density
+        ModuleIO::write_vdata_palgrid(this->Pgrid,
+                                      this->chr.rho[is],
+                                      is,
+                                      PARAM.inp.nspin,
+                                      istep,
+                                      fn,
+                                      this->pelec->eferm.get_efval(is),
+                                      &(ucell),
+                                      3,
+                                      1);
+    }
+
+    ModuleBase::timer::tick("ESolver_DM2rho", "runner");
+}
+
+template <typename TK, typename TR>
+void ESolver_DM2rho<TK, TR>::after_all_runners(UnitCell& ucell)
+{
+    ModuleBase::TITLE("ESolver_DM2rho", "after_all_runners");
+    ModuleBase::timer::tick("ESolver_DM2rho", "after_all_runners");
+
+    ESolver_KS_LCAO<TK, TR>::after_all_runners(ucell);
+
+    ModuleBase::timer::tick("ESolver_DM2rho", "after_all_runners");
+};
+
+template class ESolver_DM2rho<std::complex<double>, double>;
+template class ESolver_DM2rho<std::complex<double>, std::complex<double>>;
+
+} // namespace ModuleESolver

source/module_esolver/esolver_dm2rho.h

@@ -0,0 +1,25 @@
+#ifndef ESOLVER_DM2RHO_H
+#define ESOLVER_DM2RHO_H
+
+#include "module_esolver/esolver_ks_lcao.h"
+
+#include <memory>
+
+namespace ModuleESolver
+{
+
+template <typename TK, typename TR>
+class ESolver_DM2rho : public ESolver_KS_LCAO<TK, TR>
+{
+  public:
+    ESolver_DM2rho();
+    ~ESolver_DM2rho();
+
+    void before_all_runners(UnitCell& ucell, const Input_para& inp) override;
+
+    void after_all_runners(UnitCell& ucell) override;
+
+    void runner(UnitCell& ucell, const int istep) override;
+};
+} // namespace ModuleESolver
+#endif

source/module_esolver/esolver_fp.cpp

@@ -152,20 +152,10 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
         {
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
-                double* data = nullptr;
-                if (PARAM.inp.dm_to_rho)
-                {
-                    data = this->chr.rho[is];
-                    this->pw_rhod->real2recip(this->chr.rho[is], this->chr.rhog[is]);
-                }
-                else
-                {
-                    data = this->chr.rho_save[is];
-                    this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
-                }
+                this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
                 std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
                 ModuleIO::write_vdata_palgrid(Pgrid,
-                                              data,
+                                              this->chr.rho_save[is],
                                               is,
                                               PARAM.inp.nspin,
                                               istep,
@@ -381,19 +371,16 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
     {
         if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
         {
-            std::complex<double>** rhog_tot
-                = (PARAM.inp.dm_to_rho) ? this->chr.rhog : this->chr.rhog_save;
-            double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr.rho : this->chr.rho_save;
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
-                this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]);
+                this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
             }
             ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-CHARGE-DENSITY.restart",
                                  PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local,
                                  this->pw_rhod,
                                  PARAM.inp.nspin,
                                  ucell.GT,
-                                 rhog_tot,
+                                 this->chr.rhog_save,
                                  GlobalV::MY_POOL,
                                  GlobalV::RANK_IN_POOL,
                                  GlobalV::NPROC_IN_POOL);

source/module_esolver/esolver_ks.cpp

@@ -439,13 +439,6 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
     // 2) before_scf (electronic iteration loops)
     this->before_scf(ucell, istep);
 
-    // 3) write charge density
-    if (PARAM.inp.dm_to_rho)
-    {
-        ModuleBase::timer::tick(this->classname, "runner");
-        return; // nothing further is needed
-    }
-
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
     // 4) SCF iterations

source/module_esolver/esolver_ks_lcao.cpp

@@ -637,13 +637,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
                                          this->pelec->nelec_spin,
                                          this->pelec->skip_weights);
 
-            if (!PARAM.inp.dm_to_rho)
-            {
-                auto _pelec = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-                elecstate::calEBand(_pelec->ekb,_pelec->wg,_pelec->f_en);
-                elecstate::cal_dm_psi(_pelec->DM->get_paraV_pointer(), _pelec->wg, *this->psi, *(_pelec->DM));
-                _pelec->DM->cal_DMR();
-            }
+            auto _pelec = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
+            elecstate::calEBand(_pelec->ekb, _pelec->wg, _pelec->f_en);
+            elecstate::cal_dm_psi(_pelec->DM->get_paraV_pointer(), _pelec->wg, *this->psi, *(_pelec->DM));
+            _pelec->DM->cal_DMR();
+
             this->pelec->psiToRho(*this->psi);
             this->pelec->skip_weights = false;
 

source/module_esolver/lcao_before_scf.cpp

@@ -251,57 +251,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->cal_DMR();
     }
 
-    if (PARAM.inp.dm_to_rho)
-    {
-        // file name of DM
-        std::string zipname = "output_DM0.npz";
-        elecstate::DensityMatrix<TK, double>* dm
-            = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
-
-        // read DM from file
-        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
-
-        // if nspin=2, need extra reading
-        if (PARAM.inp.nspin == 2)
-        {
-            zipname = "output_DM1.npz";
-            ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
-        }
-
-        elecstate::calculate_weights(this->pelec->ekb,
-                                     this->pelec->wg,
-                                     this->pelec->klist,
-                                     this->pelec->eferm,
-                                     this->pelec->f_en,
-                                     this->pelec->nelec_spin,
-                                     this->pelec->skip_weights);
-
-        this->pelec->psiToRho(*this->psi);
-
-        int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
-
-        for (int is = 0; is < nspin0; is++)
-        {
-            std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
-
-            // write electron density
-            ModuleIO::write_vdata_palgrid(this->Pgrid,
-                                          this->chr.rho[is],
-                                          is,
-                                          PARAM.inp.nspin,
-                                          istep,
-                                          fn,
-                                          this->pelec->eferm.get_efval(is),
-                                          &(ucell),
-                                          3,
-                                          1);
-        }
-
-        // why we need to return here? mohan add 2025-03-10
-        ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf");
-        return;
-    }
-
     // 16) the electron charge density should be symmetrized,
     // here is the initialization
     Symmetry_rho srho;

source/module_hsolver/hsolver_lcao.cpp

@@ -83,13 +83,11 @@ void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
                                      pes->f_en,
                                      pes->nelec_spin,
                                      pes->skip_weights);
-        if (!PARAM.inp.dm_to_rho)
-        {
-            auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
-            elecstate::calEBand(_pes_lcao->ekb,_pes_lcao->wg,_pes_lcao->f_en);
-            elecstate::cal_dm_psi(_pes_lcao->DM->get_paraV_pointer(), _pes_lcao->wg, psi, *(_pes_lcao->DM));
-            _pes_lcao->DM->cal_DMR();
-        }
+
+        auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
+        elecstate::calEBand(_pes_lcao->ekb, _pes_lcao->wg, _pes_lcao->f_en);
+        elecstate::cal_dm_psi(_pes_lcao->DM->get_paraV_pointer(), _pes_lcao->wg, psi, *(_pes_lcao->DM));
+        _pes_lcao->DM->cal_DMR();
 
         if (!skip_charge)
         {