diff --git a/source/module_basis/module_nao/CMakeLists.txt b/source/module_basis/module_nao/CMakeLists.txt
index 1268400eac..2b848af517 100644
--- a/source/module_basis/module_nao/CMakeLists.txt
+++ b/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)
diff --git a/source/module_basis/module_nao/atomic_radials.cpp b/source/module_basis/module_nao/atomic_radials.cpp
index 73e0bb9945..5210770bb0 100644
--- a/source/module_basis/module_nao/atomic_radials.cpp
+++ b/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;
diff --git a/source/module_basis/module_nao/atomic_radials.h b/source/module_basis/module_nao/atomic_radials.h
index 84a4e87213..b4aac608f5 100644
--- a/source/module_basis/module_nao/atomic_radials.h
+++ b/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
diff --git a/source/module_basis/module_nao/beta_radials.cpp b/source/module_basis/module_nao/beta_radials.cpp
new file mode 100644
index 0000000000..782d4e0e2f
--- /dev/null
+++ b/source/module_basis/module_nao/beta_radials.cpp
@@ -0,0 +1,433 @@
+#include "module_basis/module_nao/beta_radials.h"
+
+#include "module_base/parallel_common.h"
+#include "module_base/tool_quit.h"
+
+void BetaRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank) {
+    /*
+     * Build a BetaRadials object from beta functions read from a pseudopotential file
+     *
+     * NOTE: only the rank == 0 process reads the file
+     *                                                                                      */
+    cleanup();
+
+    std::ifstream ifs;
+    bool is_open = false;
+
+    if (rank == 0)
+    {
+        ifs.open(file);
+        is_open = ifs.is_open();
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_bool(is_open);
+#endif
+
+    if (!is_open)
+    {
+        ModuleBase::WARNING_QUIT("BetaRadials::read", "Couldn't open pseudopotential file: " + file);
+    }
+
+    if (ptr_log)
+    {
+        (*ptr_log) << "\n\n\n\n";
+        (*ptr_log) << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
+        (*ptr_log) << " |                                                                   |" << std::endl;
+        (*ptr_log) << " |               SETUP BETA PROJECTORS                               |" << std::endl;
+        (*ptr_log) << " |                                                                   |" << std::endl;
+        (*ptr_log) << " | Projector information includes the cutoff radius, angular         |" << std::endl;
+        (*ptr_log) << " | momentum, zeta number and numerical values on a radial grid.      |" << std::endl;
+        (*ptr_log) << " |                                                                   |" << std::endl;
+        (*ptr_log) << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
+        (*ptr_log) << "\n\n\n\n";
+    }
+
+    itype_ = itype;
+
+    // identify the version of UPF file
+    // UPF 2.0.1 format starts with <UPF version="2.0.1">
+    int upf_version = 100;
+    if (rank == 0) {
+        std::string first_line;
+        std::getline(ifs, first_line);
+        if (first_line.find("2.0.1") != std::string::npos) {
+            upf_version = 201;
+        }
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_int(upf_version);
+#endif
+
+    switch (upf_version) {
+        case 100:
+            read_beta_upf100(ifs, ptr_log, rank);
+            break;
+        case 201:
+            read_beta_upf201(ifs, ptr_log, rank);
+            break;
+        default: /* not supposed to happen */;
+    }
+
+    if (rank == 0)
+    {
+        ifs.close();
+    }
+}
+
+void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) {
+    /*
+     * Read the nonlocal beta functions from the pseudopotential file (in old UPF format)
+     *                                                                                      */
+    std::stringstream ss;
+    std::string tmp;
+    std::string line;
+    int ngrid_max = 0;
+    bool is_good = false;
+
+    if (rank == 0) {
+        /*
+         * Read the header, including
+         *
+         * 1. Element symbol
+         * 2. Maximum angular momentum
+         * 3. Number of radial grid points
+         * 4. Number of beta functions
+         *                                                                                  */
+        while (std::getline(ifs, line)) {
+            if (line.find("Element") != std::string::npos) {
+                ss.str("");
+                ss << line;
+                ss >> symbol_;
+            } else if (line.find("Max angular momentum component") != std::string::npos) {
+                ss.str("");
+                ss << line;
+                ss >> lmax_;
+            } else if (line.find("Number of points in mesh") != std::string::npos) {
+                ss.str("");
+                ss << line;
+                ss >> ngrid_max;
+            } else if (line.find("Number of Projectors") != std::string::npos) {
+                ss.str("");
+                ss << line;
+                ss >> tmp >> nchi_; // nchi_ is the number of beta functions
+            } 
+
+            // should obtain valid nchi_, symbol_ & ngrid_max upon reaching the end of header
+            if (line.find("</PP_HEADER>") != std::string::npos) {
+                // lmax could be -1 when there is no beta projectors, see, e.g., H.pz-vbc.UPF
+                is_good = (nchi_ >= 0) && symbol_.size() && (ngrid_max > 0);
+                break;
+            }
+        }
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_bool(is_good);
+    Parallel_Common::bcast_string(symbol_);
+    Parallel_Common::bcast_int(lmax_);
+    Parallel_Common::bcast_int(nchi_);
+    Parallel_Common::bcast_int(ngrid_max);
+#endif
+
+    if (!is_good)
+    {
+        ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf100", "PP_HEADER error");
+    }
+
+    // In case some pseudopotential file does not have any beta function
+    if (nchi_ == 0) {
+        return;
+    }
+
+    double* rgrid = new double[ngrid_max];
+    if (rank == 0) {
+        /*
+         * Read the radial grid
+         *                                                                                  */
+        while (ifs >> tmp) {
+            if (tmp == "<PP_R>") {
+                break;
+            }
+        }
+        assert(!ifs.eof());
+
+        for (int ir = 0; ir != ngrid_max; ++ir) {
+            ifs >> rgrid[ir];
+        }
+
+        ifs >> tmp;
+        assert(tmp == "</PP_R>");
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_double(rgrid, ngrid_max);
+#endif
+
+    assert(lmax_ >= 0);
+    nzeta_ = new int[lmax_ + 1];
+    for (int l = 0; l <= lmax_; ++l) {
+        nzeta_[l] = 0;
+    }
+
+    // rbeta, l & ngrid are directly read from file
+    double* rbeta = new double[ngrid_max]; // beta function is given as r*beta(r)
+    int l = 0;
+    int ngrid = 0;
+
+    int l_last = -1;
+    int izeta = 0;
+
+    chi_ = new NumericalRadial[nchi_];
+    for (int i = 0; i != nchi_; ++i) {
+        if (rank == 0) {
+            /*
+             * Read the beta functions, including
+             *
+             * 1. Angular momentum
+             * 2. Number of actual radial grid points (should be smaller than ngrid_max)
+             * 3. Numerical values (r*beta(r)) on the radial grid
+             *
+             * and record the zeta number (izeta).
+             *                                                                              */
+            while (std::getline(ifs, line)) {
+                if (line.find("<PP_BETA>") != std::string::npos) {
+                    break;
+                }
+            }
+
+            ifs >> tmp >> l;
+            ifs >> tmp >> tmp; // skip "Beta" "L"
+            ifs >> ngrid;
+
+            if (l == l_last) {
+                izeta += 1;
+            } else {
+                izeta = 0;
+                l_last = l;
+            }
+
+            for (int ir = 0; ir != ngrid; ++ir) {
+                ifs >> rbeta[ir];
+            }
+
+            nzeta_[l] += 1;
+
+            ifs >> tmp;
+            assert(tmp == "</PP_BETA>");
+        } // rank == 0
+
+#ifdef __MPI
+        Parallel_Common::bcast_int(l);
+        Parallel_Common::bcast_int(ngrid);
+        Parallel_Common::bcast_int(izeta);
+        Parallel_Common::bcast_double(rbeta, ngrid);
+#endif
+        chi_[i].build(l, true, ngrid, rgrid, rbeta, 1, izeta, symbol_, itype_);
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_int(nzeta_, lmax_ + 1);
+#endif
+
+    indexing();
+
+    delete[] rgrid;
+    delete[] rbeta;
+}
+
+void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) {
+    /*
+     * Read the nonlocal beta functions from the pseudopotential file (in UPF 2.0.1 format)
+     *                                                                                      */
+    std::string tmp;
+    std::string line;
+    int ngrid_max = 0;
+    bool is_good = false;
+
+    if (rank == 0) {
+        /*
+         * Read the header, including
+         *
+         * 1. Element symbol
+         * 2. Maximum angular momentum
+         * 3. Number of radial grid points
+         * 4. Number of beta functions
+         *                                                                                  */
+        while (std::getline(ifs, line)) {
+            if (line.find("element=") != std::string::npos) {
+                symbol_ = trim201(line);
+            } else if (line.find("l_max=") != std::string::npos) {
+                lmax_ = std::stoi(trim201(line));
+            } else if (line.find("mesh_size=") != std::string::npos) {
+                ngrid_max = std::stoi(trim201(line));
+            } else if (line.find("number_of_proj=") != std::string::npos) {
+                nchi_ = std::stoi(trim201(line));
+            }
+
+            if (line.find("/>") != std::string::npos) {
+                is_good = (nchi_ >= 0) && (ngrid_max >= 0) && symbol_.size();
+                break;
+            }
+        }
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_bool(is_good);
+    Parallel_Common::bcast_string(symbol_);
+    Parallel_Common::bcast_int(lmax_);
+    Parallel_Common::bcast_int(nchi_);
+    Parallel_Common::bcast_int(ngrid_max);
+#endif
+
+    if (!is_good)
+    {
+        ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf201", "PP_HEADER error");
+    }
+
+    // In case some pseudopotential file does not have any beta function
+    if (nchi_ == 0) {
+        return;
+    }
+
+    double* rgrid = new double[ngrid_max];
+    if (rank == 0) {
+        /*
+         * Read the radial grid
+         *                                                                                  */
+        while (std::getline(ifs, line)) {
+            if (line.find("<PP_R") != std::string::npos) {
+                break;
+            }
+        }
+        assert(!ifs.eof());
+
+        for (int ir = 0; ir != ngrid_max; ++ir) {
+            ifs >> rgrid[ir];
+        }
+
+        ifs >> tmp;
+        assert(tmp == "</PP_R>");
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_double(rgrid, ngrid_max);
+#endif
+
+    nzeta_ = new int[lmax_ + 1];
+    for (int l = 0; l <= lmax_; ++l) {
+        nzeta_[l] = 0;
+    }
+
+    // rbeta, l & ngrid are directly read from file
+    double* rbeta = new double[ngrid_max];
+    int l = 0;
+    int ngrid = 0;
+
+    int l_last = -1;
+    int izeta = 0;
+
+    chi_ = new NumericalRadial[nchi_];
+    for (int i = 0; i != nchi_; ++i) {
+        if (rank == 0) {
+            /*
+             * Read the beta functions, including
+             *
+             * 1. Angular momentum
+             * 2. Number of actual radial grid points (should be smaller than ngrid_max)
+             * 3. Numerical values on the radial grid
+             *
+             * and record the zeta number (izeta).
+             *                                                                              */
+            while (std::getline(ifs, line)) {
+                if (line.find("<PP_BETA") != std::string::npos) {
+                    break;
+                }
+            }
+            assert(!ifs.eof());
+
+            while (std::getline(ifs, line)) {
+                if (line.find("angular_momentum") != std::string::npos) {
+                    l = std::stoi(trim201(line));
+                } else if (line.find("size") != std::string::npos) {
+                    ngrid = std::stoi(trim201(line));
+                } 
+                // neither "cutoff_radius_index" nor "cutoff_radius" is reliable!
+                // the code will read all the values first and then reverse scan to determine the grid size
+
+                if (line.find(">") != std::string::npos) {
+                    is_good &= (l >= 0) && (l <= lmax_) && (ngrid > 0) && (ngrid <= ngrid_max);
+                    break;
+                }
+            }
+
+            if (l == l_last) {
+                izeta += 1;
+            } else {
+                izeta = 0;
+                l_last = l;
+            }
+
+            for (int ir = 0; ir != ngrid; ++ir) {
+                ifs >> rbeta[ir];
+            }
+
+            // reverse scan to determine the grid size
+            for (; ngrid > 0; --ngrid) {
+                if (std::abs(rbeta[ngrid-1]) > 1e-12) {
+                    break;
+                }
+            }
+
+            nzeta_[l] += 1;
+
+            ifs >> tmp;
+            assert(tmp.find("</PP_BETA") != std::string::npos);
+            is_good &= ifs.good() && (ngrid > 0);
+        } // rank == 0
+
+#ifdef __MPI
+        Parallel_Common::bcast_bool(is_good);
+        Parallel_Common::bcast_int(l);
+        Parallel_Common::bcast_int(ngrid);
+        Parallel_Common::bcast_int(izeta);
+        Parallel_Common::bcast_double(rbeta, ngrid);
+#endif
+
+        if (!is_good)
+        {
+            ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf201", "PP_BETA error");
+        }
+
+        chi_[i].build(l, true, ngrid, rgrid, rbeta, 1, izeta, symbol_, itype_);
+    }
+
+#ifdef __MPI
+    Parallel_Common::bcast_int(nzeta_, lmax_ + 1);
+#endif
+
+    indexing();
+
+    delete[] rgrid;
+    delete[] rbeta;
+}
+
+std::string BetaRadials::trim201(std::string const& str) {
+
+    // extract the substring between quotation marks (with whitespace trimmed)
+    // str MUST contain exactly a pair of quotation marks
+
+    std::string::size_type start = str.find('"');
+    std::string::size_type end = str.find_last_of('"');
+    std::string tmp = str.substr(start+1, end-start-1);
+
+    if (tmp.length() == 0) {
+        return tmp;
+    }
+
+    start = tmp.find_first_not_of(" \t");
+    end = tmp.find_last_not_of(" \t");
+    return tmp.substr(start, end+1-start);
+}
+
diff --git a/source/module_basis/module_nao/beta_radials.h b/source/module_basis/module_nao/beta_radials.h
new file mode 100644
index 0000000000..4f11dd0bf9
--- /dev/null
+++ b/source/module_basis/module_nao/beta_radials.h
@@ -0,0 +1,54 @@
+#ifndef BETA_RADIALS_H_
+#define BETA_RADIALS_H_
+
+#include "module_basis/module_nao/radial_set.h"
+
+//! The radial part of beta functions of a single element
+/*!
+ *  This class represents the radial part of all Kleinman-Bylander beta
+ *  functions of a single element as read from a pseudopotential file.
+ *
+ *  @see RadialSet
+ *
+ *  Usage:
+ *
+ *      int element_index = 1;
+ *      std::ofstream ofs_log("/path/to/log/file");
+ *      std::string upf_file = "/path/to/pseudopotential/file";
+ *
+ *      BetaRadials O_beta;
+ *      O_beta.build(orb_file, element_index, ofs_log, GlobalV::MY_RANK);
+ *
+ *                                                                          */
+class BetaRadials: public RadialSet
+{
+  public:
+    BetaRadials() {}
+    ~BetaRadials() {}
+
+    //! Build the class from a pseudopotential file
+    void build(const std::string& file, //!< pseudopotential file name
+            const int itype = 0, //!< element index in calculation
+            std::ofstream* ptr_log = nullptr, //!< output file stream for logging
+            const int rank = 0 //!< MPI rank
+    );
+
+  private:
+
+    //! Read beta projectors from a pseudopotential file of UPF 1.0.0 format
+    void read_beta_upf100(std::ifstream& ifs,               //!< input file stream from orbital file
+                         std::ofstream* ptr_log = nullptr, //!< output file stream for logging
+                         const int rank = 0                //!< MPI rank
+    );
+
+    //! Read beta projectors from a pseudopotential file of UPF 2.0.1 format
+    void read_beta_upf201(std::ifstream& ifs,               //!< input file stream from orbital file
+                         std::ofstream* ptr_log = nullptr, //!< output file stream for logging
+                         const int rank = 0                //!< MPI rank
+    );
+
+    //! extract the substring between a pair of quotation marks (for UPF v2.0.1)
+    std::string trim201(std::string const& str);
+};
+
+#endif
diff --git a/source/module_basis/module_nao/numerical_radial.cpp b/source/module_basis/module_nao/numerical_radial.cpp
index 9e0b708ef0..87b42549f0 100644
--- a/source/module_basis/module_nao/numerical_radial.cpp
+++ b/source/module_basis/module_nao/numerical_radial.cpp
@@ -11,8 +11,8 @@
 
 #include "module_base/constants.h"
 #include "module_base/cubic_spline.h"
-#include "module_base/spherical_bessel_transformer.h"
 #include "module_base/math_integral.h"
+#include "module_base/spherical_bessel_transformer.h"
 
 using ModuleBase::PI;
 
@@ -26,45 +26,45 @@ NumericalRadial::NumericalRadial()
 
 NumericalRadial::NumericalRadial(const NumericalRadial& other)
 {
-    this->symbol_ = other.symbol_;
-    this->itype_ = other.itype_;
-    this->izeta_ = other.izeta_;
-    this->l_ = other.l_;
+    symbol_ = other.symbol_;
+    itype_ = other.itype_;
+    izeta_ = other.izeta_;
+    l_ = other.l_;
 
-    this->nr_ = other.nr_;
-    this->nk_ = other.nk_;
+    nr_ = other.nr_;
+    nk_ = other.nk_;
 
-    this->is_fft_compliant_ = other.is_fft_compliant_;
+    is_fft_compliant_ = other.is_fft_compliant_;
 
-    this->pr_ = other.pr_;
-    this->pk_ = other.pk_;
+    pr_ = other.pr_;
+    pk_ = other.pk_;
 
-    this->use_internal_transformer_ = other.use_internal_transformer_;
+    use_internal_transformer_ = other.use_internal_transformer_;
 
     // deep copy
     if (other.ptr_rgrid())
     {
-        this->rgrid_ = new double[nr_];
-        this->rvalue_ = new double[nr_];
-        std::memcpy(this->rgrid_, other.rgrid_, nr_ * sizeof(double));
-        std::memcpy(this->rvalue_, other.rvalue_, nr_ * sizeof(double));
+        rgrid_ = new double[nr_];
+        rvalue_ = new double[nr_];
+        std::memcpy(rgrid_, other.rgrid_, nr_ * sizeof(double));
+        std::memcpy(rvalue_, other.rvalue_, nr_ * sizeof(double));
     }
 
     if (other.ptr_kgrid())
     {
-        this->kgrid_ = new double[nk_];
-        this->kvalue_ = new double[nk_];
-        std::memcpy(this->kgrid_, other.kgrid_, nk_ * sizeof(double));
-        std::memcpy(this->kvalue_, other.kvalue_, nk_ * sizeof(double));
+        kgrid_ = new double[nk_];
+        kvalue_ = new double[nk_];
+        std::memcpy(kgrid_, other.kgrid_, nk_ * sizeof(double));
+        std::memcpy(kvalue_, other.kvalue_, nk_ * sizeof(double));
     }
 
     if (use_internal_transformer_)
     {
-        this->sbt_ = new ModuleBase::SphericalBesselTransformer;
+        sbt_ = new ModuleBase::SphericalBesselTransformer;
     }
     else
     {
-        this->sbt_ = other.sbt_;
+        sbt_ = other.sbt_;
     }
 }
 
@@ -116,7 +116,8 @@ NumericalRadial& NumericalRadial::operator=(const NumericalRadial& rhs)
 
     if (rhs.use_internal_transformer_)
     {
-        if (!use_internal_transformer_) {
+        if (!use_internal_transformer_)
+        {
             sbt_ = new ModuleBase::SphericalBesselTransformer;
         }
     }
@@ -422,12 +423,12 @@ void NumericalRadial::radtab(const char op,
 
     // currently only FFT-compliant grids are supported!
     // FFT-based transform requires that two NumericalRadial objects have exactly the same grid
-    assert(this->is_fft_compliant_ && ket.is_fft_compliant_);
-    assert(this->nr_ == ket.nr_);
-    assert(this->rcut() == ket.rcut());
+    assert(is_fft_compliant_ && ket.is_fft_compliant_);
+    assert(nr_ == ket.nr_);
+    assert(rcut() == ket.rcut());
 
     double* ktmp = new double[nk_];
-    std::transform(this->kvalue_, this->kvalue_ + nk_, ket.kvalue_, ktmp, std::multiplies<double>());
+    std::transform(kvalue_, kvalue_ + nk_, ket.kvalue_, ktmp, std::multiplies<double>());
 
     int op_pk = 0;
     switch (op)
@@ -445,14 +446,14 @@ void NumericalRadial::radtab(const char op,
     { // derivative of the radial table
         if (l == 0)
         { // j'_0(x) = -j_1(x)
-            sbt_->radrfft(1, nk_, this->kcut(), ktmp, table, this->pk_ + ket.pk_ + op_pk - 1);
+            sbt_->radrfft(1, nk_, kcut(), ktmp, table, pk_ + ket.pk_ + op_pk - 1);
             std::for_each(table, table + nr_, [](double& x) { x *= -1; });
         }
         else
         { // (2*l+1) * j'_l(x) = l * j_{l-1}(x) - (l+1) * j_{l+1}(x)
             double* rtmp = new double[nr_];
-            sbt_->radrfft(l + 1, nk_, this->kcut(), ktmp, table, this->pk_ + ket.pk_ + op_pk - 1);
-            sbt_->radrfft(l - 1, nk_, this->kcut(), ktmp, rtmp, this->pk_ + ket.pk_ + op_pk - 1);
+            sbt_->radrfft(l + 1, nk_, kcut(), ktmp, table, pk_ + ket.pk_ + op_pk - 1);
+            sbt_->radrfft(l - 1, nk_, kcut(), ktmp, rtmp, pk_ + ket.pk_ + op_pk - 1);
             std::transform(table, table + nr_, rtmp, table, [l](double x, double y) {
                 return (l * y - (l + 1) * x) / (2 * l + 1);
             });
@@ -461,33 +462,38 @@ void NumericalRadial::radtab(const char op,
     }
     else
     { // radial table
-        sbt_->radrfft(l, nk_, this->kcut(), ktmp, table, this->pk_ + ket.pk_ + op_pk);
+        sbt_->radrfft(l, nk_, kcut(), ktmp, table, pk_ + ket.pk_ + op_pk);
     }
 
     delete[] ktmp;
 }
 
-void NumericalRadial::normalize(bool for_r_space) {
+void NumericalRadial::normalize(bool for_r_space)
+{
     int& ngrid = for_r_space ? nr_ : nk_;
+
+    // tbu stands for "to be updated"
     double*& grid_tbu = for_r_space ? rgrid_ : kgrid_;
     double*& value_tbu = for_r_space ? rvalue_ : kvalue_;
 
-    double factor = 0.0; 
+    double factor = 0.0;
     double* integrand = new double[ngrid];
-    double* rab = new double[ngrid-1];
+    double* rab = new double[ngrid];
 
     std::adjacent_difference(grid_tbu, grid_tbu + ngrid, rab);
-    std::transform(value_tbu, value_tbu+ngrid, grid_tbu, integrand, std::multiplies<double>());
+    std::transform(value_tbu, value_tbu + ngrid, grid_tbu, integrand, std::multiplies<double>());
     std::for_each(integrand, integrand + ngrid, [](double& x) { x *= x; });
 
+    // FIXME Simpson_Integral should use only ngrid-1 rab points!
+    rab[ngrid - 1] = rab[ngrid - 2];
     ModuleBase::Integral::Simpson_Integral(ngrid, integrand, rab, factor);
-    factor = 1./std::sqrt(factor);
+    factor = 1. / std::sqrt(factor);
 
     std::for_each(value_tbu, value_tbu + ngrid, [factor](double& x) { x *= factor; });
     transform(for_r_space);
     delete[] rab;
     delete[] integrand;
-    //unit test TBD!
+    // unit test TBD!
 }
 
 void NumericalRadial::transform(const bool forward)
diff --git a/source/module_basis/module_nao/numerical_radial.h b/source/module_basis/module_nao/numerical_radial.h
index 88dcad87e0..9beb32d86f 100644
--- a/source/module_basis/module_nao/numerical_radial.h
+++ b/source/module_basis/module_nao/numerical_radial.h
@@ -85,7 +85,7 @@ class NumericalRadial
                const int p = 0,           //!< exponent of the implicit power term in input values, @see @ref group1
                const int izeta = 0,       //!< index for the multiplicity of radial functions of the same itype and l
                const std::string symbol = "", //!< chemical symbol
-               const int itype = 0        //!< index for the element in calculation
+               const int itype = 0            //!< index for the element in calculation
     );
 
     //! Sets a SphericalBesselTransformer.
@@ -153,9 +153,6 @@ class NumericalRadial
     //! Saves the data to file (what data, in what format?)
     void save(std::string file = "" /*! file name */) const;
 
-    //! Set itype for the object.
-    void set_itype(const int itype) { itype_ = itype; }
-
     //! Computes the radial table for two-center integrals.
     /*!
      *  TODO add support for non-FFT-compliant grid
@@ -211,100 +208,60 @@ class NumericalRadial
      *                                                                                  */
     ///@{
     //! gets symbol_
-    std::string const& symbol() const
-    {
-        return symbol_;
-    }
+    std::string const& symbol() const { return symbol_; }
 
     //! gets itype_
-    int itype() const
-    {
-        return itype_;
-    }
+    int itype() const { return itype_; }
 
     //! gets izeta_
-    int izeta() const
-    {
-        return izeta_;
-    }
+    int izeta() const { return izeta_; }
 
     //! gets the angular momentum
-    int l() const
-    {
-        return l_;
-    }
+    int l() const { return l_; }
 
     //! gets the number of r-space grid points
-    int nr() const
-    {
-        return nr_;
-    }
+    int nr() const { return nr_; }
 
     //! gets the number of k-space grid points
-    int nk() const
-    {
-        return nk_;
-    }
+    int nk() const { return nk_; }
 
     //! gets r-space grid cutoff distance
     double rcut() const
     {
-        return rgrid_ ? rgrid_[nr_ - 1] : 0.0;
+        assert(rgrid_);
+        return rgrid_[nr_ - 1];
     }
 
     //! gets k-space grid cutoff distance
     double kcut() const
     {
-        return kgrid_ ? kgrid_[nk_ - 1] : 0.0;
+        assert(kgrid_);
+        return kgrid_[nk_ - 1];
     }
 
     //! gets the pointer to r-space grid points
-    const double* ptr_rgrid() const
-    {
-        return rgrid_;
-    }
+    const double* ptr_rgrid() const { return rgrid_; }
 
     //! gets the pointer to k-space grid points
-    const double* ptr_kgrid() const
-    {
-        return kgrid_;
-    }
+    const double* ptr_kgrid() const { return kgrid_; }
 
     //! gets the pointer to r-space values
-    const double* ptr_rvalue() const
-    {
-        return rvalue_;
-    }
+    const double* ptr_rvalue() const { return rvalue_; }
 
     //! gets the pointer to k-space values
-    const double* ptr_kvalue() const
-    {
-        return kvalue_;
-    }
+    const double* ptr_kvalue() const { return kvalue_; }
 
     //! gets the exponent of the pre-multiplied power term in rvalues_. @see pr_
-    double pr() const
-    {
-        return pr_;
-    }
+    double pr() const { return pr_; }
 
     //! gets the exponent of the pre-multiplied power term in kvalues_. @see pk_
-    double pk() const
-    {
-        return pk_;
-    }
+    double pk() const { return pk_; }
 
     //! gets the flag for FFT-compliancy. @see is_fft_compliant_
-    bool is_fft_compliant() const
-    {
-        return is_fft_compliant_;
-    }
+    bool is_fft_compliant() const { return is_fft_compliant_; }
 
     //! gets the pointer to the SphericalBesselTransformer
-    const ModuleBase::SphericalBesselTransformer* ptr_sbt() const
-    {
-        return sbt_;
-    }
+    const ModuleBase::SphericalBesselTransformer* ptr_sbt() const { return sbt_; }
 
     ///@}
 
diff --git a/source/module_basis/module_nao/radial_collection.cpp b/source/module_basis/module_nao/radial_collection.cpp
new file mode 100644
index 0000000000..7131a90a34
--- /dev/null
+++ b/source/module_basis/module_nao/radial_collection.cpp
@@ -0,0 +1,95 @@
+#include "module_basis/module_nao/radial_collection.h"
+
+#include "module_base/spherical_bessel_transformer.h"
+#include "module_basis/module_nao/atomic_radials.h"
+#include "module_basis/module_nao/beta_radials.h"
+
+double RadialCollection::rcut_max() const
+{
+    double rmax = 0.0;
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        rmax = std::max(rmax, radset_[itype]->rcut_max());
+    }
+    return rmax;
+}
+
+RadialCollection::~RadialCollection()
+{
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        delete radset_[itype];
+    }
+    delete[] radset_;
+}
+
+void RadialCollection::cleanup()
+{
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        delete radset_[itype];
+    }
+    delete[] radset_;
+    radset_ = nullptr;
+    ntype_ = 0;
+    lmax_ = -1;
+    nchi_ = 0;
+}
+
+void RadialCollection::build(const int nfile, const std::string* const file, const char file_type)
+{
+    assert(file_type == 'o' || file_type == 'p');
+
+    cleanup();
+
+    ntype_ = nfile;
+    radset_ = new RadialSet*[ntype_];
+    switch (file_type)
+    {
+    case 'o':
+        for (int itype = 0; itype < ntype_; ++itype)
+        {
+            radset_[itype] = new AtomicRadials;
+            radset_[itype]->build(file[itype], itype);
+        }
+        break;
+    case 'p':
+        for (int itype = 0; itype < ntype_; ++itype)
+        {
+            radset_[itype] = new BetaRadials;
+            radset_[itype]->build(file[itype], itype);
+        }
+        break;
+    default:; /* not supposed to happen */
+    }
+
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        lmax_ = std::max(lmax_, radset_[itype]->lmax());
+        nchi_ += radset_[itype]->nchi();
+    }
+}
+
+void RadialCollection::set_transformer(ModuleBase::SphericalBesselTransformer* const sbt, const int update)
+{
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        radset_[itype]->set_transformer(sbt, update);
+    }
+}
+
+void RadialCollection::set_grid(const bool for_r_space, const int ngrid, const double* grid, const char mode)
+{
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        radset_[itype]->set_grid(for_r_space, ngrid, grid, mode);
+    }
+}
+
+void RadialCollection::set_uniform_grid(const bool for_r_space, const int ngrid, const double cutoff, const char mode)
+{
+    for (int itype = 0; itype < ntype_; ++itype)
+    {
+        radset_[itype]->set_uniform_grid(for_r_space, ngrid, cutoff, mode);
+    }
+}
diff --git a/source/module_basis/module_nao/radial_collection.h b/source/module_basis/module_nao/radial_collection.h
new file mode 100644
index 0000000000..db9e3be062
--- /dev/null
+++ b/source/module_basis/module_nao/radial_collection.h
@@ -0,0 +1,105 @@
+#ifndef RADIAL_COLLECTION_H_
+#define RADIAL_COLLECTION_H_
+
+#include <string>
+
+#include "module_basis/module_nao/radial_set.h"
+
+//! A class that holds the collection of all numerical radial functions of a kind.
+/*!
+ *  This class is supposed to hold the radial functions of all numerical atomic
+ *  orbitals or Kleinman-Bylander beta functions from all elements involved in
+ *  calculation. The way it stores the radial functions is type(element)-wise,
+ *  i.e. all radial functions of the same itype (element) are stored in a
+ *  RadialSet object, and the collection holds an array of such objects.
+ *                                                                              */
+class RadialCollection
+{
+  public:
+    RadialCollection(){};
+    ~RadialCollection();
+
+    void build(const int nfile, const std::string* const file, const char type = 'o');
+
+    /*! @name Getters
+     *
+     *  Get access to private members (and their properties).
+     *                                                                      */
+    //!@{
+    //! element symbol of a given type
+    const std::string& symbol(const int itype) const { return radset_[itype]->symbol(); }
+
+    //! number of RadialSet objects in the collection
+    int ntype() const { return ntype_; }
+
+    //! maximum angular momentum of all NumericalRadial objects in the collection
+    int lmax() const { return lmax_; }
+
+    //! maximum cutoff radius of all NumericalRadial objects in the collection
+    double rcut_max() const;
+
+    //! number of distinct radial functions of a given type and angular momentum
+    int nzeta(const int itype, const int l) const { return radset_[itype]->nzeta(l); }
+
+    //! maximum number of distinct radial functions of a given type among all angular momentum
+    int nzeta_max(const int itype) const { return radset_[itype]->nzeta_max(); }
+
+    //! total number of NumericalRadial objects in the collection
+    int nchi() const { return nchi_; }
+
+    //! number of NumericalRadial objects of a given type
+    int nchi(const int itype) const { return radset_[itype]->nchi(); }
+
+    //! get access to the NumericalRadial object with given type, angular momentum and zeta number
+    const NumericalRadial& operator()(const int itype, const int l, const int izeta) const
+    {
+        assert(itype >= 0 && itype < ntype_);
+        return radset_[itype]->chi(l, izeta);
+    }
+
+    //! get access to the RadialSet object with given type
+    const RadialSet& operator()(const int itype) const
+    {
+        assert(itype >= 0 && itype < ntype_);
+        return *radset_[itype];
+    }
+    //!@}
+
+    /*! @name property setters for all RadialSet objects
+     *
+     *  @see RadialSet
+     *                                                                     */
+    //!@{
+    //! Set a spherical Bessel transformers for all RadialSet objects
+    //! @see RadialSet::set_transformer
+    void set_transformer(ModuleBase::SphericalBesselTransformer* const sbt = nullptr, const int update = 0);
+
+    //! Set a common grid for all RadialSet objects
+    //! @see RadialSet::set_grid
+    void set_grid(const bool for_r_space, const int ngrid, const double* grid, const char mode = 'i');
+
+    //! Set a common uniform grid for all RadialSet objects
+    //! @see RadialSet::set_uniform_grid
+    void set_uniform_grid(const bool for_r_space, const int ngrid, const double cutoff, const char mode = 'i');
+    //!@}
+
+  private:
+    int ntype_ = 0; //!< number of RadialSet in the collection
+    int lmax_ = -1; //!< maximum angular momentum of all NumericalRadial objects in the collection
+    int nchi_ = 0;  //!< total number of NumericalRadial objects in the collection
+
+    //! array of RadialSet objects
+    /*!
+     *  Each object hold a set of NumericalRadial objects that belong to a single element.
+     *  The set could be either atomic orbitals or beta functions, in which case the
+     *  underlying objects are AtomicRadials or BetaRadials, respectively.
+     *
+     *  NOTE: AtomicRadials and BetaRadials do not necessarily have the same size as
+     *  RadialSet. Therefore, a multilevel pointer is necessary for polymorphism.
+     *                                                                                    */
+    RadialSet** radset_ = nullptr;
+
+    void cleanup();
+};
+
+#endif
diff --git a/source/module_basis/module_nao/radial_set.cpp b/source/module_basis/module_nao/radial_set.cpp
index 4dab92210a..f194370175 100644
--- a/source/module_basis/module_nao/radial_set.cpp
+++ b/source/module_basis/module_nao/radial_set.cpp
@@ -1,5 +1,7 @@
 #include "module_basis/module_nao/radial_set.h"
 
+#include <algorithm>
+
 RadialSet::~RadialSet()
 {
     delete[] nzeta_;
@@ -7,21 +9,147 @@ RadialSet::~RadialSet()
     delete[] index_map_;
 }
 
-int RadialSet::index(int l, int izeta) const
+double RadialSet::rcut_max() const
+{
+    double rmax = 0.0;
+    for (int i = 0; i != nchi_; ++i)
+    {
+        if (chi_[i].rcut() > rmax)
+        {
+            rmax = chi_[i].rcut();
+        }
+    }
+    return rmax;
+}
+
+RadialSet::RadialSet(const RadialSet& other)
+{
+    symbol_ = other.symbol_;
+    itype_ = other.itype_;
+    lmax_ = other.lmax_;
+
+    nzeta_ = nullptr;
+    if (lmax_ >= 0)
+    {
+        nzeta_ = new int[lmax_ + 1];
+        for (int l = 0; l <= lmax_; l++)
+            nzeta_[l] = other.nzeta_[l];
+    }
+    nzeta_max_ = other.nzeta_max_;
+    nchi_ = other.nchi_;
+
+    chi_ = nullptr;
+    if (nchi_ > 0)
+    {
+        chi_ = new NumericalRadial[nchi_];
+        for (int i = 0; i < nchi_; i++)
+        {
+            chi_[i] = other.chi_[i]; // deep copy
+        }
+    }
+
+    index_map_ = nullptr;
+    if (lmax_ >= 0 && nzeta_max_ > 0)
+    {
+        index_map_ = new int[(lmax_ + 1) * nzeta_max_];
+        for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++)
+        {
+            index_map_[i] = other.index_map_[i];
+        }
+    }
+}
+
+RadialSet& RadialSet::operator=(const RadialSet& rhs)
+{
+    if (&rhs == this)
+    {
+        return *this;
+    }
+
+    symbol_ = rhs.symbol_;
+    itype_ = rhs.itype_;
+    lmax_ = rhs.lmax_;
+
+    delete[] nzeta_;
+    nzeta_ = nullptr;
+    if (lmax_ >= 0)
+    {
+        nzeta_ = new int[lmax_ + 1];
+        for (int l = 0; l <= lmax_; l++)
+            nzeta_[l] = rhs.nzeta_[l];
+    }
+    nzeta_max_ = rhs.nzeta_max_;
+    nchi_ = rhs.nchi_;
+
+    delete[] chi_;
+    chi_ = nullptr;
+    if (nchi_ > 0)
+    {
+        chi_ = new NumericalRadial[nchi_];
+        for (int i = 0; i < nchi_; i++)
+        {
+            chi_[i] = rhs.chi_[i]; // deep copy
+        }
+    }
+
+    delete[] index_map_;
+    index_map_ = nullptr;
+    if (lmax_ >= 0 && nzeta_max_ > 0)
+    {
+        index_map_ = new int[(lmax_ + 1) * nzeta_max_];
+        for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++)
+        {
+            index_map_[i] = rhs.index_map_[i];
+        }
+    }
+    return *this;
+}
+
+int RadialSet::index(const int l, const int izeta) const
 {
     assert(l >= 0 && l <= lmax_);
     assert(izeta >= 0 && izeta < nzeta_[l]);
     return index_map_[l * nzeta_max_ + izeta];
 }
 
-const NumericalRadial& RadialSet::chi(int l, int izeta)
+void RadialSet::indexing()
+{
+    if (!nzeta_)
+    {
+        return;
+    }
+
+    assert(lmax_ >= 0);
+    nzeta_max_ = *std::max_element(nzeta_, nzeta_ + lmax_ + 1);
+
+    delete[] index_map_;
+    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;
+            }
+        }
+    }
+}
+
+const NumericalRadial& RadialSet::chi(const int l, const int izeta)
 {
     int i = index_map_[l * nzeta_max_ + izeta];
     assert(i >= 0 && i < nchi_);
     return chi_[i];
 }
 
-void RadialSet::set_transformer(ModuleBase::SphericalBesselTransformer* sbt, int update)
+void RadialSet::set_transformer(ModuleBase::SphericalBesselTransformer* const sbt, const int update)
 {
     for (int i = 0; i < nchi_; i++)
         chi_[i].set_transformer(sbt, update);
@@ -50,8 +178,6 @@ void RadialSet::cleanup()
     nzeta_max_ = 0;
     nchi_ = 0;
 
-    rcut_max_ = 0;
-
     delete[] chi_;
     chi_ = nullptr;
 
diff --git a/source/module_basis/module_nao/radial_set.h b/source/module_basis/module_nao/radial_set.h
index acc28e9c3e..7b1b1b1417 100644
--- a/source/module_basis/module_nao/radial_set.h
+++ b/source/module_basis/module_nao/radial_set.h
@@ -7,22 +7,24 @@
 
 #include "module_basis/module_nao/numerical_radial.h"
 
-//! An abstract class representing a related set of numerical radial functions
+//! An abstract class representing a related set of numerical radial functions.
 /*!
  *  This abstract class represents a set of numerical radial functions from
- *  a single file and of the same "itype" number. This is supposed to be the
- *  base class for concrete classes like AtomicRadials and BetaRadials, in
- *  which case "itype" is the element index in calculation, and they represent
- *  the set of all radial functions of the numerical atomic orbitals and
- *  Kleinman-Bylander beta functions of a single element respectively.
+ *  a single file. This is supposed to be the base class for concrete classes
+ *  like AtomicRadials and BetaRadials, in which case they represent the set
+ *  of all radial functions of the numerical atomic orbitals and Kleinman-
+ *  Bylander beta functions of a single element respectively.
  *
  *  @see AtomicRadials BetaRadials
  *                                                                          */
 class RadialSet
 {
   public:
-    RadialSet(){};
-    ~RadialSet();
+    RadialSet() {}
+    RadialSet(const RadialSet&);            //!< deep copy
+    RadialSet& operator=(const RadialSet&); //!< deep copy
+
+    virtual ~RadialSet();
 
     //! build the set of numerical radial functions from a file
     virtual void build(const std::string& file,                //!< orbital or pseudopotential file
@@ -37,38 +39,20 @@ class RadialSet
      *  Get access to private members.
      *                                                                      */
     //!@{
-    const std::string& symbol() const
-    {
-        return symbol_;
-    }
-    int itype() const
-    {
-        return itype_;
-    }
-    int lmax() const
-    {
-        return lmax_;
-    }
+    const std::string& symbol() const { return symbol_; }
+    int itype() const { return itype_; }
+    int lmax() const { return lmax_; }
+    double rcut_max() const;
 
-    int nzeta(int l) const
+    int nzeta(const int l) const
     {
+        assert(l >= 0 && l <= lmax_);
         return nzeta_[l];
     }
-    int nzeta_max() const
-    {
-        return nzeta_max_;
-    }
-    int nchi() const
-    {
-        return nchi_;
-    }
-
-    double rcut_max() const
-    {
-        return rcut_max_;
-    }
+    int nzeta_max() const { return nzeta_max_; }
+    int nchi() const { return nchi_; }
 
-    const NumericalRadial& chi(int l, int izeta);
+    const NumericalRadial& chi(const int l, const int izeta);
     //!@}
 
     /*! @name property setters for all NumericalRadial objects
@@ -78,7 +62,7 @@ class RadialSet
     //!@{
     //! Set a spherical Bessel transformers for all NumericalRadial objects
     //! @see NumericalRadial::set_transformer
-    void set_transformer(ModuleBase::SphericalBesselTransformer* sbt = nullptr, int update = 0);
+    void set_transformer(ModuleBase::SphericalBesselTransformer* const sbt = nullptr, const int update = 0);
 
     //! Set a common grid for all NumericalRadial objects
     //! @see NumericalRadial::set_grid
@@ -98,8 +82,6 @@ class RadialSet
     int nzeta_max_ = 0;    //!< maximum number of NumericalRadial objects among each angular momentum
     int nchi_ = 0;         //!< total number of NumericalRadial objects
 
-    double rcut_max_ = 0; //!< maximum r-space cutoff radius among all NumericalRadial objects
-
     NumericalRadial* chi_ = nullptr; //!< array of NumericalRadial objects
 
     int* index_map_ = nullptr; //!< map (l,izeta) to an index in chi_ array
@@ -109,6 +91,9 @@ class RadialSet
 
     //! get the index in chi_ array from (l,izeta)
     int index(const int l, const int izeta) const;
+
+    //! calculate nzeta_max_ and build index_map_ from nzeta_ and lmax_
+    void indexing();
 };
 
 #endif
diff --git a/source/module_basis/module_nao/test/CMakeLists.txt b/source/module_basis/module_nao/test/CMakeLists.txt
index 108abc0c8b..a530d7fe69 100644
--- a/source/module_basis/module_nao/test/CMakeLists.txt
+++ b/source/module_basis/module_nao/test/CMakeLists.txt
@@ -1,5 +1,3 @@
-remove_definitions(-D__MPI)
-
 AddTest(
   TARGET numerical_radial
   SOURCES
@@ -10,7 +8,7 @@ AddTest(
 
 AddTest(
   TARGET atomic_radials
-  SOURCES 
+  SOURCES
     atomic_radials_test.cpp
     ../atomic_radials.cpp
     ../radial_set.cpp
@@ -18,3 +16,24 @@ AddTest(
   LIBS ${math_libs} device base
 )
 
+AddTest(
+  TARGET beta_radials
+  SOURCES
+    beta_radials_test.cpp
+    ../beta_radials.cpp
+    ../radial_set.cpp
+    ../numerical_radial.cpp
+  LIBS ${math_libs} device base
+)
+
+AddTest(
+  TARGET radial_collection
+  SOURCES
+    radial_collection_test.cpp
+    ../radial_collection.cpp
+    ../atomic_radials.cpp
+    ../beta_radials.cpp
+    ../radial_set.cpp
+    ../numerical_radial.cpp
+  LIBS ${math_libs} device base
+)
diff --git a/source/module_basis/module_nao/test/atomic_radials_test.cpp b/source/module_basis/module_nao/test/atomic_radials_test.cpp
index f215b2c2a2..d493dc7ed2 100644
--- a/source/module_basis/module_nao/test/atomic_radials_test.cpp
+++ b/source/module_basis/module_nao/test/atomic_radials_test.cpp
@@ -7,6 +7,7 @@
 #endif
 
 #include "module_base/constants.h"
+#include "module_base/global_variable.h"
 
 /***********************************************************
  *      Unit test of class "AtomicRadials"
@@ -14,7 +15,7 @@
 /*!
  *  Tested functions:
  *
- *  - read
+ *  - build
  *      - parse an orbital file and initialize the NumericalRadial objects
  *
  *  - all "getters"
@@ -38,8 +39,11 @@ class AtomicRadialsTest : public ::testing::Test
 
 TEST_F(AtomicRadialsTest, ReadAndGet)
 {
+#ifdef __MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK);
+#endif
 
-    Ti_radials.build(file);
+    Ti_radials.build(file, 0, nullptr, GlobalV::MY_RANK);
 
     EXPECT_EQ(Ti_radials.lmax(), 3);
     EXPECT_EQ(Ti_radials.nzeta(0), 4);
@@ -69,9 +73,12 @@ TEST_F(AtomicRadialsTest, ReadAndGet)
 
 TEST_F(AtomicRadialsTest, BatchSet)
 {
+#ifdef __MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK);
+#endif
 
     int itype = 5;
-    Ti_radials.build(file, itype);
+    Ti_radials.build(file, itype, nullptr, GlobalV::MY_RANK);
 
     EXPECT_EQ(Ti_radials.itype(), 5);
     EXPECT_EQ(Ti_radials.chi(0, 0).itype(), 5);
diff --git a/source/module_basis/module_nao/test/beta_radials_test.cpp b/source/module_basis/module_nao/test/beta_radials_test.cpp
new file mode 100644
index 0000000000..056dacf2c5
--- /dev/null
+++ b/source/module_basis/module_nao/test/beta_radials_test.cpp
@@ -0,0 +1,195 @@
+#include "module_basis/module_nao/beta_radials.h"
+
+#include "gtest/gtest.h"
+
+#ifdef __MPI
+#include <mpi.h>
+#endif
+
+#include "module_base/constants.h"
+#include "module_base/global_variable.h"
+
+/***********************************************************
+ *      Unit test of class "BetaRadials"
+ ***********************************************************/
+/*!
+ *  Tested functions:
+ *
+ *  - build
+ *      - parse a pseudopotential file and initialize individual NumericalRadial objects
+ *
+ *  - all "getters"
+ *      - Get access to private members.
+ *
+ *  - all "batch setters"
+ *      - Set a property for all NumericalRadial objects at once
+ *                                                                      */
+class BetaRadialsTest : public ::testing::Test
+{
+  protected:
+    void SetUp();
+    void TearDown() {}
+
+    BetaRadials beta; //!< object under test
+
+    std::string dir = "../../../../../tests/PP_ORB/"; //!< directory with test files
+    std::string log_file = "./test_files/atomic_orbital.log"; //!< file for logging
+};
+
+void BetaRadialsTest::SetUp() {
+#ifdef __MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK);
+#endif
+}
+
+TEST_F(BetaRadialsTest, ReadAndGet100)
+{
+    /*
+     * Read beta projectors from a UPF file of old format
+     *                                                                  */
+    std::string file100 = "Zn.LDA.UPF"; //!< a UPF 1.0.0 file to read from
+    beta.build(dir+file100, 314, nullptr, GlobalV::MY_RANK);
+
+    EXPECT_EQ(beta.itype(), 314);
+    EXPECT_EQ(beta.symbol(), "Zn");
+    EXPECT_EQ(beta.lmax(), 2);
+    EXPECT_EQ(beta.nzeta(0), 0);
+    EXPECT_EQ(beta.nzeta(1), 1);
+    EXPECT_EQ(beta.nzeta(2), 1);
+    EXPECT_EQ(beta.nzeta_max(), 1);
+    EXPECT_EQ(beta.nchi(), 2);
+
+    EXPECT_DOUBLE_EQ(beta.rcut_max(),  3.36331676102);
+
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[0]  , 3.21829939971e-05);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[4]  , 3.39007851980e-05);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[888], 3.31987661585e+00);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[889], 3.36331676102e+00);
+
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[0]  , 3.21829939971e-05);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[4]  , 3.39007851980e-05);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[887], 3.27699753773e+00);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[888], 3.31987661585e+00);
+
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[0]  , -9.73759791529e-09);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[4]  , -1.08048430719e-08);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[600], -5.66090228695e-02);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[888], -1.89079314582e-11);
+    EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[889],  1.99483148995e-12);
+
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[0]  , -2.84316719619e-11);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[4]  , -3.32316831459e-11);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[600], -3.87224488590e-01);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[887],  3.31894955664e-11);
+    EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[888], -3.22751710952e-12);
+}
+
+TEST_F(BetaRadialsTest, ReadAndGetVoid) {
+    /*
+     * This test reads a UPF file with no beta projectors.
+     * The test is to check that the code does not crash.
+     *                                                                      */
+    std::string file0 = "H.pz-vbc.UPF"; //!< a UPF file with no beta projectors
+    beta.build(dir+file0, 5, nullptr, GlobalV::MY_RANK);
+
+    EXPECT_EQ(beta.nchi(), 0);
+    EXPECT_EQ(beta.itype(), 5);
+    EXPECT_EQ(beta.lmax(), -1);
+    EXPECT_EQ(beta.symbol(), "H");
+}
+
+TEST_F(BetaRadialsTest, ReadAndGet201)
+{
+    /*
+     * This test read beta projectors from a UPF file of 2.0.1 format
+     *                                                                      */
+    std::string file201 = "Pb_ONCV_PBE-1.0.upf"; //!< a UPF 2.0.1 file to read from
+    beta.build(dir+file201, 999, nullptr, GlobalV::MY_RANK);
+
+    EXPECT_EQ(beta.itype(), 999);
+    EXPECT_EQ(beta.symbol(), "Pb");
+    EXPECT_EQ(beta.lmax(), 3);
+    EXPECT_EQ(beta.nzeta(0), 2);
+    EXPECT_EQ(beta.nzeta(1), 2);
+    EXPECT_EQ(beta.nzeta(2), 2);
+    EXPECT_EQ(beta.nzeta(3), 2);
+    EXPECT_EQ(beta.nzeta_max(), 2);
+    EXPECT_EQ(beta.nchi(), 8);
+
+    // NOTE: neither "cutoff_radius_index" nor "cutoff_radius" is reliable!
+    // the code reads all the values first and then reverse scan to determine the grid size
+    EXPECT_DOUBLE_EQ(beta.rcut_max(), 3.68);
+
+    EXPECT_EQ(beta.chi(0,0).rcut(), 3.64);
+    EXPECT_EQ(beta.chi(0,0).nr(), 365);
+    EXPECT_EQ(beta.chi(0,0).izeta(), 0);
+    EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rgrid()[0]  , 0.0000);
+    EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rgrid()[8]  , 0.0800);
+    EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rgrid()[364], 3.6400);
+    EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rvalue()[0]  , 0.0000000000e+00);
+    EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rvalue()[4]  , 5.9689893417e-02);
+    EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rvalue()[364], -4.5888625103e-07);
+
+    EXPECT_EQ(beta.chi(3,1).rcut(), 3.68);
+    EXPECT_EQ(beta.chi(3,1).nr(), 369);
+    EXPECT_EQ(beta.chi(3,1).izeta(), 1);
+    EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rgrid()[0]  , 0.0000);
+    EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rgrid()[8]  , 0.0800);
+    EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rgrid()[368], 3.6800);
+    EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rvalue()[0]  , 0.0000000000e+00);
+    EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rvalue()[4]  , 1.7908487484e-06);
+    EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rvalue()[368], -7.0309158570e-06);
+}
+
+TEST_F(BetaRadialsTest, BatchSet)
+{
+    std::string file201 = "Pb_ONCV_PBE-1.0.upf"; //!< a UPF 2.0.1 file to read from
+    beta.build(dir+file201, 999, nullptr, GlobalV::MY_RANK);
+
+    ModuleBase::SphericalBesselTransformer sbt;
+    beta.set_transformer(&sbt);
+    for (int l = 0; l != beta.lmax(); ++l) {
+        for (int izeta = 0; izeta != beta.nzeta(l); ++izeta) {
+            EXPECT_EQ(beta.chi(l, izeta).ptr_sbt(), &sbt);
+        }
+    }
+
+    beta.set_uniform_grid(true, 2001, 20.0);
+    for (int l = 0; l != beta.lmax(); ++l) {
+        for (int izeta = 0; izeta != beta.nzeta(l); ++izeta) {
+            EXPECT_EQ(beta.chi(l, izeta).nr(), 2001);
+            EXPECT_EQ(beta.chi(l, izeta).rcut(), 20.0);
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[1500], 15.0);
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rvalue()[1500], 0.0);
+        }
+    }
+
+    double grid[5] = {0.0, 1.1, 2.2, 3.3, 4.4};
+    beta.set_grid(true, 5, grid);
+    for (int l = 0; l != beta.lmax(); ++l) {
+        for (int izeta = 0; izeta != beta.nzeta(l); ++izeta) {
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[0], 0.0);
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[1], 1.1);
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[2], 2.2);
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[3], 3.3);
+            EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[4], 4.4);
+        }
+    }
+}
+
+int main(int argc, char** argv)
+{
+
+#ifdef __MPI
+    MPI_Init(&argc, &argv);
+#endif
+
+    testing::InitGoogleTest(&argc, argv);
+    int result = RUN_ALL_TESTS();
+
+#ifdef __MPI
+    MPI_Finalize();
+#endif
+
+    return result;
+}
diff --git a/source/module_basis/module_nao/test/numerical_radial_test.cpp b/source/module_basis/module_nao/test/numerical_radial_test.cpp
index 2cbac16380..429b2bc38a 100644
--- a/source/module_basis/module_nao/test/numerical_radial_test.cpp
+++ b/source/module_basis/module_nao/test/numerical_radial_test.cpp
@@ -201,7 +201,6 @@ TEST_F(NumericalRadialTest, BuildAndGet)
     EXPECT_EQ(chi.nr(), sz);
     EXPECT_EQ(chi.nk(), 0);
     EXPECT_EQ(chi.rcut(), grid[sz - 1]);
-    EXPECT_EQ(chi.kcut(), 0);
 
     ASSERT_NE(chi.ptr_rgrid(), nullptr);
     ASSERT_NE(chi.ptr_rvalue(), nullptr);
diff --git a/source/module_basis/module_nao/test/radial_collection_test.cpp b/source/module_basis/module_nao/test/radial_collection_test.cpp
new file mode 100644
index 0000000000..c3cda85cee
--- /dev/null
+++ b/source/module_basis/module_nao/test/radial_collection_test.cpp
@@ -0,0 +1,150 @@
+#include "module_basis/module_nao/radial_collection.h"
+#include "gtest/gtest.h"
+
+#ifdef __MPI
+#include <mpi.h>
+#endif
+
+/***********************************************************
+ *      Unit test of class "RadialCollection"
+ ***********************************************************/
+/*!
+ *  Tested functions:
+ *
+ *  - build
+ *      - parse a series of orbital or pseudopotential files and
+ *        initialize AtomicRadials/BetaRadials objects accordingly.
+ *
+ *  - all "getters"
+ *      - Get access to private members.
+ *
+ *  - all "batch setters"
+ *      - Set a property for all RadialSet objects at once
+ *                                                                      */
+class RadialCollectionTest : public ::testing::Test
+{
+  protected:
+    void SetUp();
+    void TearDown();
+
+    RadialCollection orb;                                         //!< object under test
+    int nfile = 0; // number of orbital/pseudopotential files
+    std::string* file = nullptr; //!< orbitals file to read from
+    std::string log_file = "./test_files/atomic_orbital.log";         //!< file for logging
+};
+
+void RadialCollectionTest::SetUp() {
+    std::string dir = "../../../../../tests/PP_ORB/";
+    nfile = 4;
+    file = new std::string[nfile];
+    file[0] = dir + "C_gga_8au_100Ry_2s2p1d.orb";
+    file[1] = dir + "H_gga_8au_60Ry_2s1p.orb";
+    file[2] = dir + "O_gga_10au_100Ry_2s2p1d.orb";
+    file[3] = dir + "Fe_gga_9au_100Ry_4s2p2d1f.orb";
+}
+
+void RadialCollectionTest::TearDown() {
+    delete[] file;
+}
+
+TEST_F(RadialCollectionTest, BuildAndGet) {
+    orb.build(nfile, file, 'o');
+
+    EXPECT_EQ(orb.symbol(0), "C");
+    EXPECT_EQ(orb.symbol(1), "H");
+    EXPECT_EQ(orb.symbol(2), "O");
+    EXPECT_EQ(orb.symbol(3), "Fe");
+
+    EXPECT_EQ(orb.ntype(), 4);
+    EXPECT_EQ(orb.lmax(), 3);
+    EXPECT_DOUBLE_EQ(orb.rcut_max(), 10.0);
+
+    EXPECT_EQ(orb.nzeta(0,0), 2);
+    EXPECT_EQ(orb.nzeta(0,1), 2);
+    EXPECT_EQ(orb.nzeta(0,2), 1);
+
+    EXPECT_EQ(orb.nzeta(1,0), 2);
+    EXPECT_EQ(orb.nzeta(1,1), 1);
+
+    EXPECT_EQ(orb.nzeta(2,0), 2);
+    EXPECT_EQ(orb.nzeta(2,1), 2);
+    EXPECT_EQ(orb.nzeta(2,2), 1);
+
+    EXPECT_EQ(orb.nzeta(3,0), 4);
+    EXPECT_EQ(orb.nzeta(3,1), 2);
+    EXPECT_EQ(orb.nzeta(3,2), 2);
+    EXPECT_EQ(orb.nzeta(3,3), 1);
+
+    EXPECT_EQ(orb.nzeta_max(0), 2);
+    EXPECT_EQ(orb.nzeta_max(1), 2);
+    EXPECT_EQ(orb.nzeta_max(2), 2);
+    EXPECT_EQ(orb.nzeta_max(3), 4);
+
+    EXPECT_EQ(orb.nchi(0), 5);
+    EXPECT_EQ(orb.nchi(1), 3);
+    EXPECT_EQ(orb.nchi(2), 5);
+    EXPECT_EQ(orb.nchi(3), 9);
+    EXPECT_EQ(orb.nchi(), 22);
+
+    for (int itype = 0; itype <= 3; ++itype) {
+        EXPECT_EQ(orb(itype).itype(), itype);
+    }
+
+    for (int itype = 0; itype <= 3; ++itype) {
+        for (int l = 0; l <= orb(itype).lmax(); ++l) {
+            for (int izeta = 0; izeta != orb(itype).nzeta(l); ++izeta) {
+                EXPECT_EQ(orb(itype, l, izeta).l(), l);
+            }
+        }
+    }
+}
+
+TEST_F(RadialCollectionTest, BatchSet) {
+    orb.build(nfile, file, 'o');
+
+    ModuleBase::SphericalBesselTransformer sbt;
+    orb.set_transformer(&sbt);
+    orb.set_uniform_grid(true, 2001, 20.0);
+
+    EXPECT_EQ(orb.rcut_max(), 20.0);
+    for (int itype = 0; itype != orb.ntype(); ++itype) {
+        for (int l = 0; l <= orb(itype).lmax(); ++l) {
+            for (int izeta = 0; izeta != orb.nzeta(itype, l); ++izeta) {
+                EXPECT_EQ(&sbt, orb(itype, l, izeta).ptr_sbt());
+                EXPECT_DOUBLE_EQ(orb(itype, l, izeta).rcut(), 20.0);
+            }
+        }
+    }
+
+    double* grid = new double[3];
+    grid[0] = 0.0;
+    grid[1] = 1.0;
+    grid[2] = 3.14;
+
+    orb.set_grid(true, 3, grid, 'i');
+    for (int itype = 0; itype != orb.ntype(); ++itype) {
+        for (int l = 0; l <= orb(itype).lmax(); ++l) {
+            for (int izeta = 0; izeta != orb.nzeta(itype, l); ++izeta) {
+                EXPECT_DOUBLE_EQ(orb(itype, l, izeta).rcut(), 3.14);
+            }
+        }
+    }
+    delete[] grid;
+}
+
+int main(int argc, char** argv)
+{
+
+#ifdef __MPI
+    MPI_Init(&argc, &argv);
+#endif
+
+    testing::InitGoogleTest(&argc, argv);
+    int result = RUN_ALL_TESTS();
+
+#ifdef __MPI
+    MPI_Finalize();
+#endif
+
+    return result;
+}