Refactor: new class "BetaRadials" and "RadialCollection" for module_nao

source/module_basis/module_nao/CMakeLists.txt

@@ -5,6 +5,8 @@ if(ENABLE_LCAO)
 	numerical_radial.cpp
     radial_set.cpp
     atomic_radials.cpp
+    beta_radials.cpp
+    radial_collection.cpp
   )
 
   if(ENABLE_COVERAGE)

source/module_basis/module_nao/atomic_radials.cpp

@@ -9,6 +9,7 @@
 
 void AtomicRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank)
 {
+    // deallocates all arrays and reset variables
     cleanup();
 
     std::ifstream ifs;
@@ -45,6 +46,11 @@ void AtomicRadials::build(const std::string& file, const int itype, std::ofstrea
 
     itype_ = itype;
     read_abacus_orb(ifs, ptr_log, rank);
+
+    if (rank == 0)
+    {
+        ifs.close();
+    }
 }
 
 void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, const int rank)
@@ -82,10 +88,6 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log,
             {
                 ifs >> orb_ecut_;
             }
-            else if (tmp == "Cutoff(a.u.)")
-            {
-                ifs >> rcut_max_; // orbitals in an ABACUS orbital file have the same cutoff radius
-            }
             else if (tmp == "Lmax")
             {
                 ifs >> lmax_;
@@ -116,46 +118,36 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log,
          * 3. a map from (l, izeta) to 1-d array index in chi_
          *                                                                              */
         nchi_ = 0;
-        nzeta_max_ = 0;
         for (int l = 0; l <= lmax_; ++l)
         {
             nchi_ += nzeta_[l];
-            nzeta_max_ = std::max(nzeta_[l], nzeta_max_);
-        }
-
-        index_map_ = new int[(lmax_ + 1) * nzeta_max_];
-        int index_chi = 0;
-        for (int l = 0; l <= lmax_; ++l)
-        {
-            for (int izeta = 0; izeta != nzeta_max_; ++izeta)
-            {
-                if (izeta >= nzeta_[l])
-                {
-                    index_map_[l * nzeta_max_ + izeta] = -1; // -1 means no such orbital
-                }
-                else
-                {
-                    index_map_[l * nzeta_max_ + izeta] = index_chi;
-                    ++index_chi;
-                }
-            }
         }
+        indexing(); // calculate nzeta_max_ and build index_map_
     }
 
 #ifdef __MPI
     Parallel_Common::bcast_string(symbol_);
     Parallel_Common::bcast_double(orb_ecut_);
     Parallel_Common::bcast_int(lmax_);
-    Parallel_Common::bcast_int(nzeta_, lmax_ + 1);
 
     Parallel_Common::bcast_int(nchi_);
     Parallel_Common::bcast_int(nzeta_max_);
-    Parallel_Common::bcast_int(index_map_, (lmax_ + 1) * nzeta_max_);
 
     Parallel_Common::bcast_int(ngrid);
     Parallel_Common::bcast_double(dr);
 #endif
 
+    if (rank != 0)
+    {
+        nzeta_ = new int[lmax_ + 1];
+        index_map_ = new int[(lmax_ + 1) * nzeta_max_];
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_int(nzeta_, lmax_ + 1);
+    Parallel_Common::bcast_int(index_map_, (lmax_ + 1) * nzeta_max_);
+#endif
+
     double* rvalue = new double[ngrid];
     double* rgrid = new double[ngrid];
     for (int ir = 0; ir != ngrid; ++ir)
@@ -212,11 +204,6 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log,
         chi_[index(l, izeta)].normalize();
     }
 
-    if (rank == 0)
-    {
-        ifs.close();
-    }
-
     delete[] is_read;
     delete[] rvalue;
     delete[] rgrid;

source/module_basis/module_nao/atomic_radials.h

@@ -15,20 +15,16 @@
  *      int element_index = 1;
  *      std::ofstream ofs_log("/path/to/log/file");
  *      std::string orb_file = "/path/to/orbital/file";
- *      int rank = 0;
- *      MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  *
  *      AtomicRadials O_radials;
- *      O_radials.build(orb_file, element_index, ofs_log, rank);
+ *      O_radials.build(orb_file, element_index, ofs_log, GlobalV::MY_RANK);
  *
  *                                                                          */
 class AtomicRadials : public RadialSet
 {
   public:
-    AtomicRadials(){};
-    ~AtomicRadials()
-    {
-    } // ~RadialSet() is called automatically
+    AtomicRadials() {}
+    ~AtomicRadials() {} // ~RadialSet() is called automatically
 
     //! Build the class from an orbital file
     void build(const std::string& file,          //!< orbital file name
@@ -38,10 +34,7 @@ class AtomicRadials : public RadialSet
     );
 
     //! Get the energy cutoff as given by the orbital file
-    double orb_ecut() const
-    {
-        return orb_ecut_;
-    }
+    double orb_ecut() const { return orb_ecut_; }
 
   private:
     double orb_ecut_; //!< energy cutoff as given by the orbital file