feature: output H/S matrix into a sub dir in md calculation (#1258)

* feature: make dir for HS matrix output in md * refactor: add para istep for afterscf() * feature: md output HS matrix in a subdir * fix: fix typo
deepmodeling · Sep 23, 2022 · b27155b · b27155b
1 parent d2aadb3
commit b27155b
Show file tree

Hide file tree

Showing 19 changed files with 77 additions and 19 deletions.
diff --git a/source/input.cpp b/source/input.cpp
@@ -55,6 +55,7 @@ void Input::Init(const std::string &fn)
     //----------------------------------------------------------
     ModuleBase::Global_File::make_dir_out(this->suffix,
                                           this->calculation,
+                                          this->out_mat_hs2,
                                           GlobalV::MY_RANK,
                                           this->mdp.md_restart,
                                           this->out_alllog); // xiaohui add 2013-09-01

diff --git a/source/module_base/global_file.cpp b/source/module_base/global_file.cpp
@@ -19,6 +19,7 @@ namespace ModuleBase
 void ModuleBase::Global_File::make_dir_out(
     const std::string &suffix,
 	const std::string &calculation,
+    const bool &out_hs,
     const int rank,
     const bool &restart,
     const bool out_alllog)
@@ -44,6 +45,7 @@ void ModuleBase::Global_File::make_dir_out(
 
     GlobalV::global_out_dir = prefix + suffix + "/";
     GlobalV::global_stru_dir = GlobalV::global_out_dir + "STRU/";
+    GlobalV::global_matrix_dir = GlobalV::global_out_dir + "matrix_HS/";
 
 #ifdef __MPI
     MPI_Barrier(MPI_COMM_WORLD);
@@ -122,6 +124,45 @@ void ModuleBase::Global_File::make_dir_out(
 #endif
     }
 
+    // make dir for HS matrix output in md calculation
+    if(out_hs && (calculation == "md" || calculation == "sto-md"))
+    {
+        int make_dir_matrix = 0;
+        std::string command1 =  "test -d " + GlobalV::global_matrix_dir + " || mkdir " + GlobalV::global_matrix_dir;
+
+        times = 0;
+        while(times<GlobalV::NPROC)
+        {
+            if(rank==times)
+            {
+                if ( system( command1.c_str() ) == 0 )
+                {
+                    std::cout << " MAKE THE MATRIX DIR    : " << GlobalV::global_stru_dir << std::endl;
+                    make_dir_matrix = 1;
+                }
+                else
+                {
+                    std::cout << " PROC " << rank << " CAN NOT MAKE THE MATRIX DIR !!! " << std::endl;
+                    make_dir_matrix = 0;
+                }
+            }
+#ifdef __MPI
+            Parallel_Reduce::reduce_int_all(make_dir_matrix);
+#endif
+            if(make_dir_matrix>0) break;
+            ++times;
+        }
+
+#ifdef __MPI
+        if(make_dir_matrix==0)
+        {
+            std::cout << " CAN NOT MAKE THE MATRIX DIR......." << std::endl;
+            ModuleBase::QUIT();
+        }
+        MPI_Barrier(MPI_COMM_WORLD);
+#endif
+    }
+
     std::stringstream ss,ss1;
 
     // mohan add 2010-09-12

diff --git a/source/module_base/global_file.h b/source/module_base/global_file.h
@@ -23,6 +23,7 @@ namespace Global_File
 	// called in input.cpp, after reading parameters.
 	void make_dir_out(const std::string &suffix,
 		const std::string &calculation,
+        const bool &out_hs,
 		const int rank,
         const bool &restart,
 		const bool out_alllog = false); 

diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp
@@ -135,6 +135,7 @@ std::string global_epm_pseudo_card;
 std::string global_out_dir;
 std::string global_readin_dir; // zhengdy modified
 std::string global_stru_dir;
+std::string global_matrix_dir;
 
 std::ofstream ofs_running;
 std::ofstream ofs_warning;

diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h
@@ -158,6 +158,7 @@ extern std::string global_out_dir;
 extern std::string global_orbital_dir; // liuyu add 2021-08-14
 extern std::string global_readin_dir; // zhengdy modified
 extern std::string global_stru_dir;   // liuyu add 2022-05-24 for MD STRU
+extern std::string global_matrix_dir; // liuyu add 2022-09-19 for HS matrix output 
 
 extern std::ofstream ofs_running;
 extern std::ofstream ofs_warning;

diff --git a/source/module_base/tool_title.cpp b/source/module_base/tool_title.cpp
@@ -13,7 +13,7 @@ namespace ModuleBase
 //==========================================================
 void TITLE(const std::string &class_function_name)
 {
-	return;//no output
+    return;//no output
 
 #ifdef __NORMAL
     std::cout<<" ==> "<<class_function_name<<std::endl;
@@ -27,7 +27,7 @@ void TITLE(const std::string &class_function_name)
 
 void TITLE(const std::string &class_name,const std::string &function_name)
 {
-	return;//no output
+    return;//no output
 #ifdef __NORMAL
     std::cout<<" ==> "<<class_name<<"::"<<function_name<<std::endl;
 #else
@@ -41,7 +41,7 @@ void TITLE(const std::string &class_name,const std::string &function_name)
 
 void TITLE(std::ofstream &ofs,const std::string &class_name,const std::string &function_name)
 {
-	return;// no output
+    return;// no output
 #ifdef __NORMAL
     std::cout<<"\n\n ==> "<<class_name<<"::"<<function_name<<std::endl;
 #else

diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
@@ -259,7 +259,7 @@ namespace ModuleESolver
                     if(stop) break;
                 }
             }
-            afterscf();
+            afterscf(istep);
 
             ModuleBase::timer::tick(this->classname, "Run");
         }       

diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h
@@ -47,7 +47,7 @@ namespace ModuleESolver
         // Something to do after hamilt2density function in each iter loop.
         virtual void eachiterfinish(const int iter) {};
         // Something to do after SCF iterations when SCF is converged or comes to the max iter step.
-        virtual void afterscf() {};
+        virtual void afterscf(const int istep) {};
         // <Temporary> It should be replaced by a function in Hamilt Class
         virtual void updatepot(const int istep, const int iter) {};
         // choose strategy when charge density convergence achieved

diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -647,7 +647,7 @@ void ESolver_KS_LCAO::eachiterfinish(int iter)
     GlobalC::en.calculate_etot();
 }
 
-void ESolver_KS_LCAO::afterscf()
+void ESolver_KS_LCAO::afterscf(const int istep)
 {
     for (int ik = 0; ik < this->pelec->ekb.nr; ++ik)
     {
@@ -952,7 +952,7 @@ void ESolver_KS_LCAO::afterscf()
 
     if (hsolver::HSolverLCAO::out_mat_hsR)
     {
-        this->output_HS_R(); // LiuXh add 2019-07-15
+        this->output_HS_R(istep); // LiuXh add 2019-07-15
     }
 
     if(!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS)

diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h
@@ -32,7 +32,7 @@ namespace ModuleESolver
         virtual void hamilt2density(const int istep, const int iter, const double ethr) override;
         virtual void updatepot(const int istep, const int iter) override;
         virtual void eachiterfinish(const int iter) override;
-        virtual void afterscf() override;
+        virtual void afterscf(const int istep) override;
         virtual bool do_after_converge(int& iter) override;
         int two_level_step = 0;
 
@@ -56,6 +56,7 @@ namespace ModuleESolver
 
         // output subfuncs: implemented in src_io/write_HS_R.cpp
         void output_HS_R(
+            const int &istep,
             const std::string& SR_filename = "data-SR-sparse_SPIN0.csr",
             const std::string& HR_filename_up = "data-HR-sparse_SPIN0.csr",
             const std::string HR_filename_down = "data-HR-sparse_SPIN1.csr",

diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -432,7 +432,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
     }
 }
 
-void ESolver_KS_LCAO_TDDFT::afterscf()
+void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
 {
     for (int ik = 0; ik < this->pelec_td->ekb.nr; ++ik)
     {
@@ -527,7 +527,7 @@ void ESolver_KS_LCAO_TDDFT::afterscf()
 
     if (Pdiag_Double::out_mat_hsR)
     {
-        this->output_HS_R(); // LiuXh add 2019-07-15
+        this->output_HS_R(istep); // LiuXh add 2019-07-15
     }
 }
 

diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h
@@ -28,7 +28,7 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO
     virtual void hamilt2density(const int istep, const int iter, const double ethr) override;
     virtual void eachiterinit(const int istep, const int iter) override;
     virtual void updatepot(const int istep, const int iter) override;
-    virtual void afterscf() override;
+    virtual void afterscf(const int istep) override;
     void cal_edm_tddft();
 };
 

diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -432,7 +432,7 @@ namespace ModuleESolver
     }
 
 
-    void ESolver_KS_PW::afterscf()
+    void ESolver_KS_PW::afterscf(const int istep)
     {
         for(int ik=0; ik<this->pelec->ekb.nr; ++ik)
         {

diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h
@@ -32,7 +32,7 @@ namespace ModuleESolver
         virtual void eachiterinit(const int istep, const int iter) override;
         virtual void updatepot(const int istep, const int iter) override;
         virtual void eachiterfinish(const int iter) override;
-        virtual void afterscf() override;
+        virtual void afterscf(const int istep) override;
         virtual void othercalculation(const int istep)override;
 
         //temporary, this will be removed in the future;

diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp
@@ -75,7 +75,7 @@ void ESolver_SDFT_PW::eachiterfinish(int iter)
 	//print_eigenvalue(GlobalV::ofs_running);
     GlobalC::en.calculate_etot();
 }
-void ESolver_SDFT_PW::afterscf()
+void ESolver_SDFT_PW::afterscf(const int istep)
 {
     for(int ik=0; ik<this->pelec->ekb.nr; ++ik)
     {

diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h
@@ -24,7 +24,7 @@ class ESolver_SDFT_PW: public ESolver_KS_PW
     // virtual void eachiterinit(int iter) override; 
     virtual void hamilt2density(const int istep, const int iter, const double ethr) override;
     virtual void eachiterfinish(const int iter) override; 
-    virtual void afterscf() override;
+    virtual void afterscf(const int istep) override;
     virtual void postprocess() override;
 
 public:

diff --git a/source/src_io/write_HS.cpp b/source/src_io/write_HS.cpp
@@ -864,6 +864,7 @@ void HS_Matrix::save_HSR_tr(const int current_spin, LCAO_Matrix &lm)
 }
 
 void HS_Matrix::save_HSR_sparse(
+    const int &istep,
     LCAO_Matrix &lm,
     const double& sparse_threshold,
     const bool &binary,  
@@ -980,9 +981,18 @@ void HS_Matrix::save_HSR_sparse(
 
     std::stringstream ssh[2];
     std::stringstream sss;
-    ssh[0] << GlobalV::global_out_dir << HR_filename_up;
-    ssh[1] << GlobalV::global_out_dir << HR_filename_down;
-    sss << GlobalV::global_out_dir << SR_filename;
+    if(GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "sto-md")
+    {
+        ssh[0] << GlobalV::global_matrix_dir << istep << "_" << HR_filename_up;
+        ssh[1] << GlobalV::global_matrix_dir << istep << "_" << HR_filename_down;
+        sss << GlobalV::global_matrix_dir << istep << "_" << SR_filename;
+    }
+    else
+    {
+        ssh[0] << GlobalV::global_out_dir << HR_filename_up;
+        ssh[1] << GlobalV::global_out_dir << HR_filename_down;
+        sss << GlobalV::global_out_dir << SR_filename;
+    }
     std::ofstream g1[2];
     std::ofstream g2;
 

diff --git a/source/src_io/write_HS.h b/source/src_io/write_HS.h
@@ -20,6 +20,7 @@ namespace HS_Matrix
 
     // jingan add 2021-6-4, modify 2021-12-2
     void save_HSR_sparse(
+        const int &istep,
         LCAO_Matrix &lm,
         const double& sparse_threshold,
         const bool &binary,  

diff --git a/source/src_io/write_HS_R.cpp b/source/src_io/write_HS_R.cpp
@@ -10,6 +10,7 @@ namespace ModuleESolver
 // The 'sparse_threshold' is the accuracy of the sparse matrix. 
 // If the absolute value of the matrix element is less than or equal to the 'sparse_threshold', it will be ignored.
 void ESolver_KS_LCAO::output_HS_R(
+    const int &istep,
     const std::string &SR_filename,
     const std::string &HR_filename_up,
     const std::string HR_filename_down,
@@ -64,7 +65,7 @@ void ESolver_KS_LCAO::output_HS_R(
         }
     }
 
-    HS_Matrix::save_HSR_sparse(*this->UHM.LM, sparse_threshold, binary, SR_filename, HR_filename_up, HR_filename_down);
+    HS_Matrix::save_HSR_sparse(istep, *this->UHM.LM, sparse_threshold, binary, SR_filename, HR_filename_up, HR_filename_down);
     this->UHM.destroy_all_HSR_sparse();
 
     if(!GlobalV::GAMMA_ONLY_LOCAL) //LiuXh 20181011