refactor: output sparse DMR by Output_HContainer instead of output_dm…

…1.h/cpp (#4374)

* refactor: output sparse DMR by Output_HContainer instead of output_dm1.h/cpp

* fix Makefile.Objects

* [pre-commit.ci lite] apply automatic fixes

* correct refdata-DMR-sparse_SPIN0.csr in 201_NO_15_f_pseudopots

* fix bug; and remove paraV in output_hcontainer

* fix bug in UT

* avoid using pointer, and add an parameter to avoid using GlobalV::NLOCAL

* add annoation

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

source/Makefile.Objects

@@ -490,8 +490,7 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_HS_R.o\
       write_dm.o\
       output_dm.o\
-      output_dm1.o\
-      write_dm_sparse.o\
+      write_dmr.o\
       sparse_matrix.o\
       output_mulliken.o\
       output_sk.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -44,6 +44,7 @@
 //---------------------------------------------------
 
 #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
+#include "module_io/write_dmr.h"
 #include "module_io/write_wfc_nao.h"
 
 namespace ModuleESolver
@@ -253,13 +254,9 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(Input& inp, UnitCell& ucell)
     // 12) initialize the potential
     if (this->pelec->pot == nullptr)
     {
-        this->pelec->pot = new elecstate::Potential(this->pw_rhod,
-                                                    this->pw_rho,
-                                                    &GlobalC::ucell,
-                                                    &(GlobalC::ppcell.vloc),
-                                                    &(this->sf),
-                                                    &(this->pelec->f_en.etxc),
-                                                    &(this->pelec->f_en.vtxc));
+        this->pelec->pot
+            = new elecstate::Potential(this->pw_rhod, this->pw_rho, &GlobalC::ucell, &(GlobalC::ppcell.vloc),
+                                       &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc));
     }
 
 #ifdef __DEEPKS
@@ -324,13 +321,9 @@ void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
         // Initialize the potential.
         if (this->pelec->pot == nullptr)
         {
-            this->pelec->pot = new elecstate::Potential(this->pw_rhod,
-                                                        this->pw_rho,
-                                                        &GlobalC::ucell,
-                                                        &(GlobalC::ppcell.vloc),
-                                                        &(this->sf),
-                                                        &(this->pelec->f_en.etxc),
-                                                        &(this->pelec->f_en.vtxc));
+            this->pelec->pot
+                = new elecstate::Potential(this->pw_rhod, this->pw_rho, &GlobalC::ucell, &(GlobalC::ppcell.vloc),
+                                           &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc));
         }
     }
 
@@ -564,15 +557,8 @@ void ESolver_KS_LCAO<TK, TR>::init_basis_lcao(ORB_control& orb_con, Input& inp,
 #endif
 
 #ifndef USE_NEW_TWO_CENTER
-    this->orb_con.set_orb_tables(GlobalV::ofs_running,
-                                 *uot_,
-                                 GlobalC::ORB,
-                                 ucell.lat0,
-                                 GlobalV::deepks_setorb,
-                                 Lmax,
-                                 ucell.infoNL.nprojmax,
-                                 ucell.infoNL.nproj,
-                                 ucell.infoNL.Beta);
+    this->orb_con.set_orb_tables(GlobalV::ofs_running, *uot_, GlobalC::ORB, ucell.lat0, GlobalV::deepks_setorb, Lmax,
+                                 ucell.infoNL.nprojmax, ucell.infoNL.nproj, ucell.infoNL.Beta);
 #else
     two_center_bundle->tabulate();
 #endif
@@ -872,10 +858,8 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
     }
 
     // 11) compute magnetization, only for spin==2
-    GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
-                                                this->pelec->charge->nxyz,
-                                                this->pelec->charge->rho,
-                                                this->pelec->nelec_spin.data());
+    GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, this->pelec->charge->nxyz,
+                                                this->pelec->charge->rho, this->pelec->nelec_spin.data());
 
     // 12) calculate delta energy
     this->pelec->f_en.deband = this->pelec->cal_delta_eband();
@@ -1097,20 +1081,16 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     // 1) write charge difference into files for charge extrapolation
     if (GlobalV::CALCULATION != "scf")
     {
-        this->CE.save_files(istep,
-                            GlobalC::ucell,
+        this->CE.save_files(istep, GlobalC::ucell,
 #ifdef __MPI
                             this->pw_big,
 #endif
-                            this->pelec->charge,
-                            &this->sf);
+                            this->pelec->charge, &this->sf);
     }
 
     // 2) write density matrix for sparse matrix
-    if (this->LOC.out_dm1 == 1)
-    {
-        this->create_Output_DM1(istep).write();
-    }
+    ModuleIO::write_dmr(dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMR_vector(),
+                        GlobalV::NLOCAL, INPUT.out_dm1, false, GlobalV::out_app_flag, istep);
 
     // 3) write charge density
     if (GlobalV::out_chg)
@@ -1362,13 +1342,6 @@ ModuleIO::Output_DM ESolver_KS_LCAO<TK, TR>::create_Output_DM(int is, int iter)
 //! the 17th function of ESolver_KS_LCAO: create_Output_DM1
 //! mohan add 2024-05-11
 //------------------------------------------------------------------------------
-template <typename TK, typename TR>
-ModuleIO::Output_DM1 ESolver_KS_LCAO<TK, TR>::create_Output_DM1(int istep)
-{
-    const elecstate::DensityMatrix<complex<double>, double>* DM
-        = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM();
-    return ModuleIO::Output_DM1(GlobalV::NSPIN, istep, this->LOC, this->RA, this->kv, DM);
-}
 
 //------------------------------------------------------------------------------
 //! the 18th function of ESolver_KS_LCAO: create_Output_Mat_Sparse

source/module_esolver/esolver_ks_lcao.h

@@ -14,7 +14,6 @@
 #endif
 #include "module_basis/module_nao/two_center_bundle.h"
 #include "module_io/output_dm.h"
-#include "module_io/output_dm1.h"
 #include "module_io/output_mat_sparse.h"
 
 #include <memory>
@@ -116,9 +115,6 @@ class ESolver_KS_LCAO : public ESolver_KS<TK>
     /// @brief create ModuleIO::Output_DM object to output density matrix
     ModuleIO::Output_DM create_Output_DM(int is, int iter);
 
-    /// @brief create ModuleIO::Output_DM1 object to output sparse density matrix
-    ModuleIO::Output_DM1 create_Output_DM1(int istep);
-
     /// @brief create ModuleIO::Output_Mat_Sparse object to output sparse density matrix of H, S, T, r
     ModuleIO::Output_Mat_Sparse<TK> create_Output_Mat_Sparse(int istep);
 

source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp

@@ -1,9 +1,9 @@
 #include "output_hcontainer.h"
 
-#include <fstream>
-
 #include "module_io/sparse_matrix.h"
 
+#include <fstream>
+
 namespace hamilt
 {
 
@@ -12,18 +12,9 @@ namespace hamilt
  * @attention ofs should be open outside of this interface
  */
 template <typename T>
-Output_HContainer<T>::Output_HContainer(hamilt::HContainer<T>* hcontainer,
-                                        const Parallel_Orbitals* ParaV,
-                                        const UnitCell& ucell,
-                                        std::ostream& ofs,
-                                        double sparse_threshold,
-                                        int precision)
-    : _hcontainer(hcontainer),
-      _ParaV(ParaV),
-      _ucell(ucell),
-      _ofs(ofs),
-      _sparse_threshold(sparse_threshold),
-      _precision(precision)
+Output_HContainer<T>::Output_HContainer(hamilt::HContainer<T>* hcontainer, const int nrow, const int ncol,
+                                        std::ostream& ofs, double sparse_threshold, int precision)
+    : _hcontainer(hcontainer), _nrow(nrow), _ncol(ncol), _ofs(ofs), _sparse_threshold(sparse_threshold), _precision(precision)
 {
 }
 
@@ -65,7 +56,7 @@ template <typename T>
 void Output_HContainer<T>::write_single_R(int rx, int ry, int rz)
 {
     this->_hcontainer->fix_R(rx, ry, rz);
-    ModuleIO::SparseMatrix<T> sparse_matrix = ModuleIO::SparseMatrix<T>(this->_ParaV->nrow, this->_ParaV->ncol);
+    ModuleIO::SparseMatrix<T> sparse_matrix = ModuleIO::SparseMatrix<T>(_nrow, _ncol);
     sparse_matrix.setSparseThreshold(this->_sparse_threshold);
     for (int iap = 0; iap < this->_hcontainer->size_atom_pairs(); ++iap)
     {

source/module_hamilt_lcao/module_hcontainer/output_hcontainer.h

@@ -1,7 +1,6 @@
 #ifndef OUTPUT_HCONTAINER_H
 #define OUTPUT_HCONTAINER_H
 
-#include "module_cell/unitcell.h"
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
 
 namespace hamilt
@@ -14,12 +13,8 @@ template <typename T>
 class Output_HContainer
 {
   public:
-    Output_HContainer(hamilt::HContainer<T>* hcontainer,
-                      const Parallel_Orbitals* ParaV,
-                      const UnitCell& ucell,
-                      std::ostream& ofs,
-                      double sparse_threshold,
-                      int precision);
+    Output_HContainer(hamilt::HContainer<T>* hcontainer, const int nrow, const int ncol, std::ostream& ofs,
+                      double sparse_threshold, int precision);
     // write the matrices of all R vectors to the output stream
     void write();
 
@@ -37,13 +32,13 @@ class Output_HContainer
 
   private:
     hamilt::HContainer<T>* _hcontainer;
-    const UnitCell& _ucell;
-    const Parallel_Orbitals* _ParaV;
+    int _nrow=0; 
+    int _ncol=0;
     std::ostream& _ofs;
     double _sparse_threshold;
     int _precision;
 };
 
-} // namespace ModuleIO
+} // namespace hamilt
 
 #endif // OUTPUT_HCONTAINER_H

source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_output.cpp

@@ -1,9 +1,9 @@
+#include "../hcontainer.h"
 #include "../output_hcontainer.h"
+#include "module_cell/unitcell.h"
 
-#include "../hcontainer.h"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
-#include "module_cell/unitcell.h"
 
 /************************************************
  *  unit test of output_hcontainer.cpp
@@ -131,7 +131,7 @@ TEST_F(OutputHContainerTest, Write)
         HR.unfix_R();
     }
     double sparse_threshold = 0.1;
-    hamilt::Output_HContainer<double> output_HR(&HR, &ParaV, ucell, std::cout, sparse_threshold, 2);
+    hamilt::Output_HContainer<double> output_HR(&HR, ParaV.nrow, ParaV.ncol, std::cout, sparse_threshold, 2);
     // the first R
     testing::internal::CaptureStdout();
     output_HR.write(0, 1, 1);

source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_readCSR.cpp

@@ -1,12 +1,12 @@
-#include <fstream>
-
 #include "../hcontainer.h"
 #include "../output_hcontainer.h"
-#include "gmock/gmock.h"
-#include "gtest/gtest.h"
 #include "module_io/csr_reader.h"
 #include "prepare_unitcell.h"
 
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+#include <fstream>
+
 // mock functions
 #ifdef __LCAO
 InfoNonlocal::InfoNonlocal()
@@ -124,7 +124,7 @@ TEST_F(ReadHContainerTest, ReadAndOutputHContainer)
     // output SR
     std::ofstream ofs_out("SR.out");
     double sparse_threshold = 1e-10;
-    hamilt::Output_HContainer<double> output_SR(&SR, &paraV, *ucell, ofs_out, sparse_threshold, 8);
+    hamilt::Output_HContainer<double> output_SR(&SR, paraV.nrow, paraV.ncol, ofs_out, sparse_threshold, 8);
     // std::cout << SR.size_R_loop() << std::endl;
     // output_SR.write(-2, -1, 0);
     ofs_out << "STEP: " << 0 << std::endl;

source/module_io/CMakeLists.txt

@@ -57,9 +57,9 @@ if(ENABLE_LCAO)
       read_wfc_lcao.cpp
       write_wfc_nao.cpp
       write_dm.cpp
+      write_dmr.cpp
       dos_nao.cpp
       output_dm.cpp
-      output_dm1.cpp
       sparse_matrix.cpp
       file_reader.cpp
       csr_reader.cpp
@@ -75,7 +75,6 @@ if(ENABLE_LCAO)
       write_HS_R.cpp
       write_HS_sparse.cpp
       single_R_io.cpp
-      write_dm_sparse.cpp
       cal_r_overlap_R.cpp
       output_mat_sparse.cpp
   )