Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop

 into develop

docs/advanced/acceleration/cuda.md

@@ -29,7 +29,9 @@ To compile and use ABACUS in CUDA mode, you currently need to have an NVIDIA GPU
 
 Check the [Advanced Installation Options](https://abacus-rtd.readthedocs.io/en/latest/advanced/install.html#build-with-cuda-support) for the installation of CUDA version support.
 
-When the compilation parameter USE_ELPA is ON (which is the default value) and USE_CUDA is also set to ON, the ELPA library needs to [enable GPU support](https://github.com/marekandreas/elpa/blob/master/documentation/INSTALL.md) at compile time.
+Setting both USE_ELPA and USE_CUDA to ON does not automatically enable ELPA to run on GPUs. ELPA support for GPUs needs to be enabled when ELPA is compiled. [enable GPU support](https://github.com/marekandreas/elpa/blob/master/documentation/INSTALL.md).
+
+The ABACUS program will automatically determine whether the current ELPA supports GPU based on the elpa/elpa_configured_options.h header file. Users can also check this header file to determine the GPU support of ELPA in their environment. ELPA introduced a new API elpa_setup_gpu in version 2023.11.001. So if you want to enable ELPA GPU in ABACUS, the ELPA version must be greater than or equal to 2023.11.001.
 
 ## Run with the GPU support by editing the INPUT script:
 

docs/advanced/input_files/input-main.md

@@ -666,6 +666,7 @@ These variables are used to control parameters related to input files.
 - **Type**: String
 - **Description**: the name of the structure file
   - Containing various information about atom species, including pseudopotential files, local orbitals files, cell information, atom positions, and whether atoms should be allowed to move.
+  - When [calculation](#calculation) is set to `md` and [md_restart](#md_restart) is set to `true`, this keyword will NOT work.
   - Refer to [Doc](https://github.com/deepmodeling/abacus-develop/blob/develop/docs/advanced/input_files/stru.md)
 - **Default**: STRU
 

source/Makefile.Objects

@@ -523,6 +523,7 @@ OBJS_IO=input_conv.o\
     read_input_item_other.o\
     read_input_item_output.o\
     read_set_globalv.o\
+    orb_io.o\
 
 OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_orb_info.o\

source/driver_run.cpp

@@ -40,7 +40,7 @@ void Driver::driver_run() {
 
     // the life of ucell should begin here, mohan 2024-05-12
     // delete ucell as a GlobalC in near future
-    GlobalC::ucell.setup_cell(PARAM.inp.stru_file, GlobalV::ofs_running);
+    GlobalC::ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
     Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell,
                                          PARAM.inp.min_dist_coef);
 

source/module_base/global_variable.cpp

@@ -55,8 +55,6 @@ int GSIZE = DSIZE;
 //----------------------------------------------------------
 // EXPLAIN : The input file name and directory
 //----------------------------------------------------------
-std::string stru_file = "STRU";
-
 std::ofstream ofs_running;
 std::ofstream ofs_warning;
 std::ofstream ofs_info;   // output math lib info

source/module_base/global_variable.h

@@ -80,7 +80,6 @@ extern int KPAR_LCAO;
 // NAME : ofs_running( contain information during runnnig)
 // NAME : ofs_warning( contain warning information, including error)
 //==========================================================
-extern std::string stru_file;
 // extern std::string global_pseudo_type; // mohan add 2013-05-20 (xiaohui add
 // 2013-06-23)
 extern std::ofstream ofs_running;

source/module_base/sph_bessel_recursive.h

@@ -7,7 +7,7 @@
 #define SPH_BESSEL_RECURSIVE_H
 
 #include <vector>
-#include"stddef.h"
+#include<cstddef>
 
 namespace ModuleBase
 {
@@ -32,7 +32,7 @@ class Sph_Bessel_Recursive::D1
 
 private:
 	std::vector<std::vector<double>> jlx;		// jlx[l][x]
-	double dx;
+	double dx = 0.0;
 	bool finish_set_dx = false;
 
 	void cal_jlx_0( const int l_size );
@@ -56,7 +56,7 @@ class Sph_Bessel_Recursive::D2
 
 private:
 	std::vector<std::vector<std::vector<double>>> jlx;		// jlx[l][x1][x2]
-	double dx;
+	double dx = 0.0;
 	bool finish_set_dx = false;
 
 	void cal_jlx_0( const int l_size, const size_t ix1_size, const size_t ix2_size );

source/module_basis/module_nao/atomic_radials.cpp

@@ -3,11 +3,16 @@
 #include "module_base/math_integral.h"
 #include "module_base/parallel_common.h"
 #include "module_base/tool_quit.h"
+
+// FIXME: should update with pyabacus
+// #include "module_io/orb_io.h"
+
 #include "projgen.h"
 
 #include <fstream>
 #include <iostream>
 #include <string>
+#include <numeric>
 
 AtomicRadials& AtomicRadials::operator=(const AtomicRadials& rhs)
 {
@@ -207,7 +212,6 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log,
     {
         rgrid[ir] = ir * dr;
     }
-
     chi_ = new NumericalRadial[nchi_];
 
     // record whether an orbital has been read or not
@@ -269,116 +273,47 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log,
     delete[] rgrid;
 }
 
-void AtomicRadials::read_abacus_orb(std::ifstream& ifs,
-                                    std::string& elem,
-                                    double& ecut,
-                                    int& nr,
-                                    double& dr,
-                                    std::vector<int>& nzeta,
-                                    std::vector<std::vector<double>>& radials,
-                                    const int rank)
-{
-    nr = 0; // number of grid points
-    dr = 0; // grid spacing
-    int lmax, nchi = 0; // number of radial functions
-    std::vector<std::vector<int>> radial_map_; // build a map from [l][izeta] to 1-d array index
-    std::string tmp;
-    // first read the header
-    if (rank == 0)
-    {
-        if (!ifs.is_open())
-        {
-            ModuleBase::WARNING_QUIT("AtomicRadials::read_abacus_orb", "Couldn't open orbital file.");
-        }
-        while (ifs >> tmp)
-        {
-            if (tmp == "Element")
-            {
-                ifs >> elem;
-            }
-            else if (tmp == "Cutoff(Ry)")
-            {
-                ifs >> ecut;
-            }
-            else if (tmp == "Lmax")
-            {
-                ifs >> lmax;
-                nzeta.resize(lmax + 1);
-                for (int l = 0; l <= lmax; ++l)
-                {
-                    ifs >> tmp >> tmp >> tmp >> nzeta[l];
-                }
-            }
-            else if (tmp == "Mesh")
-            {
-                ifs >> nr;
-                continue;
-            }
-            else if (tmp == "dr")
-            {
-                ifs >> dr;
-                break;
-            }
-        }
-        radial_map_.resize(lmax + 1);
-        for (int l = 0; l <= lmax; ++l)
-        {
-            radial_map_[l].resize(nzeta[l]);
-        }
-        int ichi = 0;
-        for (int l = 0; l <= lmax; ++l)
-        {
-            for (int iz = 0; iz < nzeta[l]; ++iz)
-            {
-                radial_map_[l][iz] = ichi++; // return the value of ichi, then increment
-            }
-        }
-        nchi = ichi; // total number of radial functions
-        radials.resize(nchi);
-        std::for_each(radials.begin(), radials.end(), [nr](std::vector<double>& v) { v.resize(nr); });
-    }
-
-    // broadcast the header information
-#ifdef __MPI
-    Parallel_Common::bcast_string(elem);
-    Parallel_Common::bcast_double(ecut);
-    Parallel_Common::bcast_int(lmax);
-    Parallel_Common::bcast_int(nchi);
-    Parallel_Common::bcast_int(nr);
-    Parallel_Common::bcast_double(dr);
-#endif
-
-    // then adjust the size of the vectors
-    if (rank != 0)
-    {
-        nzeta.resize(lmax + 1);
-        radials.resize(nchi);
-        std::for_each(radials.begin(), radials.end(), [nr](std::vector<double>& v) { v.resize(nr); });
-    }
-    // broadcast the number of zeta functions for each angular momentum
-#ifdef __MPI
-    Parallel_Common::bcast_int(nzeta.data(), lmax + 1);
-#endif
-
-    // read the radial functions by rank0
-    int ichi = 0;
-    for (int i = 0; i != nchi; ++i)
-    {
-        if (rank == 0)
-        {
-            int l, izeta;
-            ifs >> tmp >> tmp >> tmp;
-            ifs >> tmp >> l >> izeta;
-            ichi = radial_map_[l][izeta];
-            for (int ir = 0; ir != nr; ++ir)
-            {
-                ifs >> radials[ichi][ir];
-            }
-        }
-    // broadcast the radial functions
-#ifdef __MPI
-        Parallel_Common::bcast_int(ichi); // let other ranks know where to store the radial function
-        Parallel_Common::bcast_double(radials[ichi].data(), nr);
-#endif
-    }
-}
+// FIXME: should update with pyabacus
+// void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, const int rank)
+// {
+//     /*
+//      * Read the orbital file.
+//      *
+//      * For orbital file format, see
+//      * (new) abacus-develop/tools/SIAB/PyTorchGradient/source/IO/print_orbital.py
+//      * (old) abacus-develop/tools/SIAB/SimulatedAnnealing/source/src_spillage/Plot_Psi.cpp
+//      *                                                                                  */
+//     int ngrid = 0; // number of grid points
+//     double dr = 0; // grid spacing
+//     std::string tmp;
+
+//     int nr;
+//     std::vector<int> nzeta;
+//     std::vector<std::vector<double>> radials;
+
+//     ModuleIO::read_abacus_orb(ifs, symbol_, orb_ecut_, nr, dr, nzeta, radials, rank);
+
+//     lmax_ = nzeta.size() - 1;
+//     nzeta_ = new int[lmax_ + 1];
+//     std::copy(nzeta.begin(), nzeta.end(), nzeta_);
+
+//     nchi_ = std::accumulate(nzeta.begin(), nzeta.end(), 0);
+//     nzeta_max_ = *std::max_element(nzeta.begin(), nzeta.end());
+
+//     indexing();
+
+//     std::vector<double> rgrid(nr);
+//     std::iota(rgrid.begin(), rgrid.end(), 0);
+//     std::for_each(rgrid.begin(), rgrid.end(), [dr](double& r) { r *= dr; });
+//     chi_ = new NumericalRadial[nchi_];
+//     int ichi = 0;
+//     for (int l = 0; l <= lmax_; ++l)
+//     {
+//         for (int izeta = 0; izeta < nzeta[l]; ++izeta)
+//         {
+//             chi_[index(l, izeta)].build(l, true, nr, rgrid.data(), radials[ichi].data(), 0, izeta, symbol_, itype_, false);
+//             chi_[index(l, izeta)].normalize();
+//             ++ichi;
+//         }
+//     }
+// }

source/module_basis/module_nao/atomic_radials.h

@@ -43,28 +43,6 @@ class AtomicRadials : public RadialSet
     //! Get the energy cutoff as given by the orbital file
     double orb_ecut() const { return orb_ecut_; }
 
-    /**
-     * @brief static version of read_abacus_orb. A delete-new operation may cause the memory leak, 
-     * it is better to use std::vector to replace the raw pointer.
-     * 
-     * @param ifs [in] ifstream from the orbital file, via `std::ifstream ifs(forb);`
-     * @param elem [out] element symbol
-     * @param ecut [out] planewave energy cutoff
-     * @param nr [out] number of radial grid points
-     * @param dr [out] radial grid spacing
-     * @param nzeta [out] number of zeta functions for each angular momentum
-     * @param radials [out] radial orbitals
-     * @param rank [in] MPI rank
-     */
-    static void read_abacus_orb(std::ifstream& ifs,
-                                std::string& elem,
-                                double& ecut,
-                                int& nr,
-                                double& dr,
-                                std::vector<int>& nzeta,
-                                std::vector<std::vector<double>>& radials,
-                                const int rank = 0);
-
   private:
     double orb_ecut_; //!< energy cutoff as given by the orbital file
 

source/module_basis/module_nao/radial_set.cpp

@@ -7,6 +7,9 @@
 
 #include "module_base/spherical_bessel_transformer.h"
 
+// FIXME: should update with pyabacus
+// #include "module_io/orb_io.h"
+
 RadialSet::~RadialSet()
 {
     delete[] nzeta_;
@@ -252,3 +255,28 @@ void RadialSet::write_abacus_orb(const std::string& file_name, const int rank) c
     }
     file_to.close();
 }
+
+// FIXME: should update with pyabacus
+// void RadialSet::write_abacus_orb(const std::string& forb, const int rank) const
+// {
+//     std::ofstream ofs;
+//     ofs.open(forb, std::ios::out);
+
+//     const double dr = 0.01;
+//     const int nr = static_cast<int>(rcut_max_ / dr) + 1;
+//     std::vector<int> nzeta(lmax_ + 1);
+//     std::copy(nzeta_, nzeta_ + lmax_ + 1, nzeta.begin());
+//     std::vector<std::vector<double>> radials(nchi_, std::vector<double>(nr, 0.0));
+//     int ichi = 0;
+//     for (int l = 0; l <= lmax_; l++)
+//     {
+//         for (int izeta = 0; izeta < nzeta[l]; izeta++)
+//         {
+//             std::copy(chi_[index(l, izeta)].rvalue(), chi_[index(l, izeta)].rvalue() + nr, radials[ichi].begin());
+//             ichi++;
+//         }
+//     }
+
+//     ModuleIO::write_abacus_orb(ofs, symbol_, 100, nr, dr, nzeta, radials, rank);
+//     ofs.close();
+// }