Feature: output R matrix at each MD step (#1825)

* Use template to reconstruct parse_expression

* Feature: output R matrix at each MD step

* Modify'matrix_HS' to 'matrix' for R matrix output

---------

Co-authored-by: dyzheng <[email protected]>

Author: 1041176461

source/module_base/global_file.cpp

@@ -20,6 +20,7 @@ void ModuleBase::Global_File::make_dir_out(
     const std::string &suffix,
 	const std::string &calculation,
     const bool &out_hs,
+    const bool &out_r,
     const int rank,
     const bool &restart,
     const bool out_alllog)
@@ -43,7 +44,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/";
+    GlobalV::global_matrix_dir = GlobalV::global_out_dir + "matrix/";
 
 #ifdef __MPI
     MPI_Barrier(MPI_COMM_WORLD);
@@ -123,7 +124,7 @@ void ModuleBase::Global_File::make_dir_out(
     }
 
     // make dir for HS matrix output in md calculation
-    if(out_hs && calculation == "md")
+    if((out_hs || out_r) && calculation == "md")
     {
         int make_dir_matrix = 0;
         std::string command1 =  "test -d " + GlobalV::global_matrix_dir + " || mkdir " + GlobalV::global_matrix_dir;

source/module_base/global_file.h

@@ -24,6 +24,7 @@ namespace Global_File
 	void make_dir_out(const std::string &suffix,
 		const std::string &calculation,
         const bool &out_hs,
+		const bool &out_r,
 		const int rank,
         const bool &restart,
 		const bool out_alllog = false); 

source/module_base/global_variable.h

@@ -163,7 +163,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::string global_matrix_dir; // liuyu add 2022-09-19 for HS matrix outpu, jiyy modified 2023-01-23 for R matrix output
 
 extern std::ofstream ofs_running;
 extern std::ofstream ofs_warning;

source/module_esolver/esolver_ks_lcao.cpp

@@ -906,11 +906,11 @@ void ESolver_KS_LCAO::afterscf(const int istep)
 
         if (hsolver::HSolverLCAO::out_mat_hsR)
         {
-            r_matrix.out_rR_other(this->LM.output_R_coor);
+            r_matrix.out_rR_other(istep, this->LM.output_R_coor);
         }
         else
         {
-            r_matrix.out_rR();
+            r_matrix.out_rR(istep);
         }
     }
 

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -423,11 +423,11 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
 
         if (hsolver::HSolverLCAO::out_mat_hsR)
         {
-            r_matrix.out_rR_other(this->LM.output_R_coor);
+            r_matrix.out_rR_other(istep, this->LM.output_R_coor);
         }
         else
         {
-            r_matrix.out_rR();
+            r_matrix.out_rR(istep);
         }
     }
 }

source/module_io/cal_r_overlap_R.cpp

@@ -226,7 +226,7 @@ void cal_r_overlap_R::init(const Parallel_Orbitals &pv)
 }
 
 
-void cal_r_overlap_R::out_rR()
+void cal_r_overlap_R::out_rR(const int &istep)
 {	
     ModuleBase::TITLE("cal_r_overlap_R", "out_rR");
     ModuleBase::timer::tick("cal_r_overlap_R", "out_rR");
@@ -406,7 +406,10 @@ void cal_r_overlap_R::out_rR()
     {
         std::ofstream out_r;
         std::stringstream ssr;
-        ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
+        if(GlobalV::CALCULATION == "md")
+            ssr << GlobalV::global_matrix_dir << istep << "_" << "data-rR-sparse.csr";
+        else
+            ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
 
         if (binary)
         {
@@ -443,7 +446,7 @@ void cal_r_overlap_R::out_rR()
 }
 
 
-void cal_r_overlap_R::out_rR_other(const std::set<Abfs::Vector3_Order<int>> &output_R_coor)
+void cal_r_overlap_R::out_rR_other(const int &istep, const std::set<Abfs::Vector3_Order<int>> &output_R_coor)
 {	
     ModuleBase::TITLE("cal_r_overlap_R", "out_rR_other");
     ModuleBase::timer::tick("cal_r_overlap_R", "out_rR_other");
@@ -456,7 +459,10 @@ void cal_r_overlap_R::out_rR_other(const std::set<Abfs::Vector3_Order<int>> &out
 
     std::ofstream out_r;
     std::stringstream ssr;
-    ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
+    if(GlobalV::CALCULATION == "md")
+        ssr << GlobalV::global_matrix_dir << istep << "_" << "data-rR-sparse.csr";
+    else
+        ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
 
     if (GlobalV::DRANK == 0)
     {

source/module_io/cal_r_overlap_R.h

@@ -34,8 +34,8 @@ class cal_r_overlap_R
     bool binary = false;
 
     void init(const Parallel_Orbitals &pv);
-    void out_rR();
-    void out_rR_other(const std::set<Abfs::Vector3_Order<int>> &output_R_coor);
+    void out_rR(const int &istep);
+    void out_rR_other(const int &istep, const std::set<Abfs::Vector3_Order<int>> &output_R_coor);
 
 private:
     void initialize_orb_table();

source/module_io/input.cpp

@@ -57,6 +57,7 @@ void Input::Init(const std::string &fn)
     ModuleBase::Global_File::make_dir_out(this->suffix,
                                           this->calculation,
                                           this->out_mat_hs2,
+                                          this->out_mat_r,
                                           GlobalV::MY_RANK,
                                           this->mdp.md_restart,
                                           this->out_alllog); // xiaohui add 2013-09-01

source/module_io/input_conv.cpp

@@ -27,7 +27,8 @@
 #include "module_elecstate/potentials/gatefield.h"
 #include "module_psi/kernels/device.h"
 
-void Input_Conv::parse_expression(const std::string &fn, std::vector<double> &vec)
+template <typename T>
+void Input_Conv::parse_expression(const std::string &fn, std::vector<T> &vec)
 {
     ModuleBase::TITLE("Input_Conv", "parse_expression");
     ModuleBase::timer::tick("Input_Conv", "parse_expression");
@@ -73,7 +74,7 @@ void Input_Conv::parse_expression(const std::string &fn, std::vector<double> &ve
         {
             int pos = sub_str.find("*");
             int num = stoi(sub_str.substr(0, pos));
-            double occ = stof(sub_str.substr(pos + 1, sub_str.size()));
+            T occ = stof(sub_str.substr(pos + 1, sub_str.size()));
             // std::vector<double> ocp_temp(num, occ);
             // const std::vector<double>::iterator dest = vec.begin() + count;
             // copy(ocp_temp.begin(), ocp_temp.end(), dest);
@@ -85,7 +86,11 @@ void Input_Conv::parse_expression(const std::string &fn, std::vector<double> &ve
         {
             // vec[count] = stof(sub_str);
             // count += 1;
-            vec.emplace_back(stof(sub_str));
+            std::stringstream convert;
+            convert << sub_str;
+            T occ;
+            convert >> occ;
+            vec.emplace_back(occ);
         }
     }
 }

source/module_io/input_conv.h

@@ -24,7 +24,8 @@ namespace Input_Conv
     // Function `parse_expression` is used to parse input parameters as expressions into vectors
     // fn (string): expressions such as "3*1 0 2*0.5 3*0"
     // arr (vector): stores parsing results, for example, "3*1 0 2*0.5 1*1.5" can be parsed as [1, 1, 1, 0, 0.5, 0.5, 1.5]
-    void parse_expression(const std::string &fn, std::vector<double> &arr);
+    template<typename T>
+    void parse_expression(const std::string &fn, std::vector<T> &arr);
 }