Refactor ModuleIO::read/write_cube() (#5150)

* Refactor ModuleIO::read/write_cube()

* Refactor ModuleIO::write_cube_core()

* Refactor ModuleIO::read_cube_core()

* add annotation in ModuleIO::read/write_cube_core()

source/module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp

@@ -182,7 +182,7 @@ void Parallel_Grid::z_distribution(void)
 
 
 #ifdef __MPI
-void Parallel_Grid::zpiece_to_all(double *zpiece, const int &iz, double *rho)
+void Parallel_Grid::zpiece_to_all(double *zpiece, const int &iz, double *rho) const
 {
 	if(PARAM.inp.esolver_type == "sdft")
 	{
@@ -256,7 +256,7 @@ void Parallel_Grid::zpiece_to_all(double *zpiece, const int &iz, double *rho)
 #endif
 
 #ifdef __MPI
-void Parallel_Grid::zpiece_to_stogroup(double *zpiece, const int &iz, double *rho)
+void Parallel_Grid::zpiece_to_stogroup(double *zpiece, const int &iz, double *rho) const
 {
 	assert(allocate);	
 	//TITLE("Parallel_Grid","zpiece_to_all");

source/module_hamilt_pw/hamilt_pwdft/parallel_grid.h

@@ -21,8 +21,8 @@ class Parallel_Grid
 		const int &nczp, const int &nrxx, const int &nbz, const int &bz); //LiuXh add 20180606
 
 #ifdef __MPI	
-	void zpiece_to_all(double *zpiece, const int &iz, double *rho);
-	void zpiece_to_stogroup(double *zpiece, const int &iz, double *rho); //qainrui add for sto-dft 2021-7-21
+	void zpiece_to_all(double *zpiece, const int &iz, double *rho) const;
+	void zpiece_to_stogroup(double *zpiece, const int &iz, double *rho) const; //qainrui add for sto-dft 2021-7-21
 
 	void reduce_to_fullrho(double *rhotot, double *rhoin);
 #endif

source/module_io/cube_io.h

@@ -10,43 +10,89 @@ namespace ModuleIO
 {
 bool read_cube(
 #ifdef __MPI
-    Parallel_Grid* Pgrid,
+    const Parallel_Grid*const Pgrid,
 #endif
-    int my_rank,
-    std::string esolver_type,
-    int rank_in_stogroup,
-    const int& is,
+    const int my_rank,
+    const std::string esolver_type,
+    const int rank_in_stogroup,
+    const int is,
     std::ofstream& ofs_running,
-    const int& nspin,
+    const int nspin,
     const std::string& fn,
-    double* data,
-    const int& nx,
-    const int& ny,
-    const int& nz,
+    double*const data,
+    const int nx,
+    const int ny,
+    const int nz,
     double& ef,
-    const UnitCell* ucell,
+    const UnitCell*const ucell,
     int& prenspin,
-    const bool& warning_flag = true);
+    const bool warning_flag = true);
 
 void write_cube(
 #ifdef __MPI
-    const int& bz,
-    const int& nbz,
-    const int& nplane,
-    const int& startz_current,
+    const int bz,
+    const int nbz,
+    const int nplane,
+    const int startz_current,
 #endif
-    const double* data,
-    const int& is,
-    const int& nspin,
-    const int& iter,
+    const double*const data,
+    const int is,
+    const int nspin,
+    const int iter,
     const std::string& fn,
-    const int& nx,
-    const int& ny,
-    const int& nz,
-    const double& ef,
-    const UnitCell* ucell,
-    const int& precision = 11,
-    const int& out_fermi = 1); // mohan add 2007-10-17
+    const int nx,
+    const int ny,
+    const int nz,
+    const double ef,
+    const UnitCell*const ucell,
+    const int precision = 11,
+    const int out_fermi = 1); // mohan add 2007-10-17
+
+
+// when MPI:
+//      read file as order (ixy,iz) to data[ixy*nz+iz]
+// when serial:
+//      read file as order (ixy,iz) to data[iz*nxy+ixy]
+void read_cube_core_match(
+    std::ifstream &ifs,
+#ifdef __MPI
+    const Parallel_Grid*const Pgrid,
+    const bool flag_read_rank,
+#endif
+    double*const data,
+    const int nxy,
+    const int nz);
+
+void read_cube_core_mismatch(
+    std::ifstream &ifs,
+#ifdef __MPI
+    const Parallel_Grid*const Pgrid,
+    const bool flag_read_rank,
+#endif
+    double*const data,
+    const int nx,
+    const int ny,
+    const int nz,
+    const int nx_read,
+    const int ny_read,
+    const int nz_read);
+
+// when MPI:
+//      write data[ixy*nplane+iz] to file as order (ixy,iz)
+// when serial:
+//      write data[iz*nxy+ixy] to file as order (ixy,iz)
+void write_cube_core(
+    std::ofstream &ofs_cube,
+#ifdef __MPI
+    const int bz,
+    const int nbz,
+    const int nplane,
+    const int startz_current,
+#endif
+    const double*const data,
+    const int nxy,
+    const int nz,
+    const int n_data_newline);
 
     /**
      * @brief The trilinear interpolation method
@@ -86,7 +132,7 @@ void write_cube(
                                const int& ny,
                                const int& nz,
 #ifdef __MPI
-                               double** data
+                               std::vector<std::vector<double>> &data
 #else
                                double* data
 #endif