Merge branch 'develop' into para19

.github/workflows/test.yml

@@ -39,7 +39,7 @@ jobs:
             --from-ref ${{ github.event.pull_request.base.sha }}
               --to-ref ${{ github.event.pull_request.head.sha }}
         continue-on-error: true
-      - uses: pre-commit-ci/[email protected].2
+      - uses: pre-commit-ci/[email protected].3
 
       - name: Build
         run: |

source/module_base/CMakeLists.txt

@@ -73,6 +73,7 @@ if(BUILD_TESTING)
     add_subdirectory(kernels/test)
     add_subdirectory(module_mixing/test)
     add_subdirectory(module_device/test)
+    add_subdirectory(grid/test)
     if (USE_ABACUS_LIBM)
       add_subdirectory(libm/test)
     endif()

source/module_base/global_variable.cpp

@@ -63,15 +63,6 @@ std::ofstream ofs_info;   // output math lib info
 std::ofstream ofs_device; // output device info
 
 
-// added by zhengdy-soc
-bool NONCOLIN = false;
-bool LSPINORB = false;
-bool DOMAG = false;
-bool DOMAG_Z = false;
-int NPOL = 1;
-
-std::vector<std::string> rpa_orbitals;
-
 //==========================================================
 // device flags added by denghui
 //==========================================================

source/module_base/global_variable.h

@@ -28,12 +28,6 @@ extern bool use_uspp;
 extern std::string KS_SOLVER;  // xiaohui add 2013-09-01
 extern double SEARCH_RADIUS;   // 11.1 // mohan add 2011-03-10
 
-// added by zhengdy-soc
-extern bool NONCOLIN;     // 0 : collinear ; 1 : non-collinear
-extern bool LSPINORB;     // 0 : no soc ; 1 : has soc
-extern bool DOMAG;        // 1 : calculate the magnetism with x, y, z component
-extern bool DOMAG_Z;      // 1 : constrain the magnetism to z axis
-extern int NPOL;          // 1 : no soc; 2 : has soc
 
 extern int NB2D;           // 16.5 dividsion of 2D_matrix.
 
@@ -94,8 +88,6 @@ extern std::ofstream ofs_warning;
 extern std::ofstream ofs_info;
 extern std::ofstream ofs_device;
 
-// rpa related
-extern std::vector<std::string> rpa_orbitals;
 
 // mixing parameters
 

source/module_base/grid/delley.h

@@ -0,0 +1,57 @@
+#ifndef GRID_DELLEY_H
+#define GRID_DELLEY_H
+
+#include <vector>
+#include <functional>
+
+/**
+ * @brief Delley's grid for quadrature on the unit sphere.
+ *
+ * Reference:
+ * Delley, B. (1996). High order integration schemes on the unit sphere.
+ * Journal of computational chemistry, 17(9), 1152-1155.
+ *
+ */
+namespace Grid {
+namespace Delley {
+
+/**
+ * @brief Number of Delley's grid points for a certain order of accuracy.
+ *
+ * This function finds the minimum Delley's grid with an accuracy
+ * of order "lmax". On exit, lmax is set to the order of this grid, and
+ * its corresponding number of grid points is returned. If no such grid
+ * is available, lmax is left unchanged and the function will return -1.
+ *
+ * For example, if lmax = 20 on input, the function will return 194 and
+ * lmax will be set to 23.
+ *
+ */
+int ngrid(int& lmax);
+
+
+/**
+ * @brief Delley's quadrature grid and weights.
+ *
+ * This function retrieves the minimum Delley's grid with an accuray
+ * of order "lmax". On exit, lmax is set to the order of this grid, and
+ * the coordinates & weights are returned in "grid" & "weight".
+ *
+ * Coordinates are stored in the following order:
+ *
+ *      x0, y0, z0, x1, y1, z1, x2, y2, z2, ...
+ *
+ * "grid" and "weight" must be pre-allocated to hold 3*ngrid(lmax) and
+ * ngrid(lmax) elements, respectively. The function will return 0 if
+ * successful, or -1 if the requested order cannot be fulfilled.
+ *
+ */
+int get(int& lmax, double* grid, double* weight);
+
+// a handy wrapper doing the same as above
+int get(int& lmax, std::vector<double>& grid, std::vector<double>& weight);
+
+} // end of namespace Delley
+} // end of namespace Grid
+
+#endif

source/module_base/grid/test/CMakeLists.txt

@@ -0,0 +1,9 @@
+remove_definitions(-D__MPI)
+
+AddTest(
+  TARGET test_delley
+  SOURCES test_delley.cpp
+  ../delley.cpp
+  ../../ylm.cpp
+)
+

source/module_base/grid/test/test_delley.cpp

@@ -0,0 +1,138 @@
+#include "module_base/grid/delley.h"
+#include "module_base/ylm.h"
+
+#include "gtest/gtest.h"
+#include <random>
+#ifdef __MPI
+#include <mpi.h>
+#endif
+
+using namespace Grid;
+
+// mock the function to prevent unnecessary dependency
+namespace ModuleBase {
+void WARNING_QUIT(const std::string&, const std::string&) {}
+}
+
+class DelleyTest: public ::testing::Test {
+protected:
+    void randgen(int lmax, std::vector<double>& coef);
+    const double tol = 1e-12;
+};
+
+
+void DelleyTest::randgen(int lmax, std::vector<double>& coef) {
+    coef.resize((lmax + 1) * (lmax + 1));
+
+    // fill coef with uniformly distributed random numbers
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<double> dis(0.0, 1.0);
+    for (size_t i = 0; i < coef.size(); ++i) {
+        coef[i] = dis(gen);
+    }
+
+    // normalize the coefficients
+    double fac = 0.0;
+    for (size_t i = 0; i < coef.size(); ++i) {
+        fac += coef[i] * coef[i];
+    }
+
+    fac = 1.0 / std::sqrt(fac);
+    for (size_t i = 0; i < coef.size(); ++i) {
+        coef[i] *= fac;
+    }
+}
+
+
+TEST_F(DelleyTest, NumGrid) {
+    int lmax = 5;
+    int ngrid = Delley::ngrid(lmax);
+    EXPECT_EQ(lmax, 17);
+    EXPECT_EQ(ngrid, 110);
+
+    lmax = 17;
+    ngrid = Delley::ngrid(lmax);
+    EXPECT_EQ(lmax, 17);
+    EXPECT_EQ(ngrid, 110);
+
+    lmax = 20;
+    ngrid = Delley::ngrid(lmax);
+    EXPECT_EQ(lmax, 23);
+    EXPECT_EQ(ngrid, 194);
+
+    lmax = 59;
+    ngrid = Delley::ngrid(lmax);
+    EXPECT_EQ(lmax, 59);
+    EXPECT_EQ(ngrid, 1202);
+
+    lmax = 60;
+    ngrid = Delley::ngrid(lmax);
+    EXPECT_EQ(lmax, 60);
+    EXPECT_EQ(ngrid, -1);
+}
+
+
+TEST_F(DelleyTest, Accuracy) {
+    /* 
+     * Given
+     *
+     *      f = c[0]*Y00 + c[1]*Y10 + c[2]*Y11 + ...,
+     *
+     * where c[0], c[1], c[2], ... are some random numbers, the integration
+     * of |f|^2 on the unit sphere
+     *
+     *      \int |f|^2 d\Omega = c[0]^2 + c[1]^2 + c[2]^2 + ... .
+     *
+     * This test verifies with the above integral that quadrature with
+     * Delley's grid is exact up to floating point errors.
+     *
+     */
+    std::vector<double> grid, weight, coef;
+
+    for (int grid_lmax = 17; grid_lmax < 60; grid_lmax +=6) {
+        Delley::get(grid_lmax, grid, weight);
+        int func_lmax = grid_lmax / 2;
+        randgen(func_lmax, coef);
+
+        double val = 0.0;
+        std::vector<double> ylm_real;
+        for (size_t i = 0; i < weight.size(); i++) {
+            ModuleBase::Ylm::sph_harm(func_lmax,
+                    grid[3*i], grid[3*i+1], grid[3*i+2], ylm_real);
+            double tmp = 0.0;
+            for (size_t i = 0; i < coef.size(); ++i) {
+                tmp += coef[i] * ylm_real[i];
+            }
+            val += weight[i] * tmp * tmp;
+        }
+        val *= 4.0 * std::acos(-1.0);
+
+        double val_ref = 0.0;
+        for (size_t i = 0; i < coef.size(); ++i) {
+            val_ref += coef[i] * coef[i];
+        }
+
+        double abs_diff = std::abs(val - val_ref);
+        EXPECT_LT(abs_diff, tol);
+        //printf("order = %2i    val_ref = %8.5f    abs_diff = %8.5e\n",
+        //        grid_lmax, val_ref, abs_diff);
+    }
+}
+
+
+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;
+}

source/module_base/module_device/device.h

@@ -44,7 +44,7 @@ int get_device_kpar(const int& kpar);
 
 /**
  * @brief Get the device flag object
- * for module_io GlobalV::device_flag
+ * for module_io PARAM.globalv.device_flag
  */
 std::string get_device_flag(const std::string& device,
                             const std::string& ks_solver,

source/module_base/ylm.cpp

@@ -54,7 +54,7 @@ std::vector<double> Ylm::ylmcoef = {
 // here Lmax == max angular momentum + 1
 void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3<double> &vec, double ylmr[] )
 {
-	ModuleBase::timer::tick ("Ylm","get_ylm_real");
+	//ModuleBase::timer::tick ("Ylm","get_ylm_real");
 	//1e-9 is too large
 	const double cut0 = 1e-12;
 	// allocate space.
@@ -165,7 +165,7 @@ void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3<double> &vec,
 		}
 	}// end do
 
-	ModuleBase::timer::tick ("Ylm", "get_ylm_real");
+	//ModuleBase::timer::tick ("Ylm", "get_ylm_real");
 	return;
 }
 
@@ -536,7 +536,7 @@ void Ylm::sph_harm
 	std::vector<double> &rly
 )
 {
-	rly.reserve( (Lmax+1)*(Lmax+1) );
+	rly.resize( (Lmax+1)*(Lmax+1) );
 
 	//begin calculation
 	/***************************

source/module_basis/module_nao/beta_radials.cpp

@@ -1,5 +1,6 @@
 #include "module_basis/module_nao/beta_radials.h"
 
+#include "module_parameter/parameter.h"
 #include "module_base/global_variable.h"
 #include "module_base/parallel_common.h"
 #include "module_base/tool_quit.h"
@@ -402,7 +403,7 @@ void BetaRadials::build(const Numerical_Nonlocal& nl, const int itype, std::ofst
 //#endif
 //
 //    // It is an error if lspinorb is set to true but the pseudopotential file does not contain spin-orbit information
-//    if (!has_so && GlobalV::LSPINORB)
+//    if (!has_so && PARAM.inp.lspinorb)
 //    {
 //        ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf201",
 //                "lspinorb is set to true but the pseudopotential file does not contain spin-orbit information");
@@ -558,7 +559,7 @@ void BetaRadials::build(const Numerical_Nonlocal& nl, const int itype, std::ofst
 //    double* rbeta_final = nullptr;
 //    if (rank == 0)
 //    {
-//        if (!GlobalV::LSPINORB && has_so)
+//        if (!PARAM.inp.lspinorb && has_so)
 //        {
 //            /*
 //             * read the PP_DIJ block for averaging
@@ -727,7 +728,7 @@ void BetaRadials::build(const Numerical_Nonlocal& nl, const int itype, std::ofst
 //    delete[] ngrid_final;
 //    delete[] rbeta_final;
 //
-//    if (rank == 0 && !GlobalV::LSPINORB && has_so)
+//    if (rank == 0 && !PARAM.inp.lspinorb && has_so)
 //    {
 //        delete[] l;
 //        delete[] ngrid;

source/module_cell/atom_pseudo.cpp

@@ -1,6 +1,7 @@
 #include "atom_pseudo.h"
 #include "module_parameter/parameter.h"
 
+#include "module_parameter/parameter.h"
 Atom_pseudo::Atom_pseudo()
 {
 }
@@ -68,7 +69,7 @@ void Atom_pseudo::set_d_so(ModuleBase::ComplexMatrix& d_so_in,
 
         if (this->lmax > -1)
         {
-            if (GlobalV::LSPINORB)
+            if (PARAM.inp.lspinorb)
             {
                 int is = 0;
                 for (int is1 = 0; is1 < 2; is1++)

source/module_cell/read_atoms.cpp

@@ -139,18 +139,6 @@ int UnitCell::read_atom_species(std::ifstream &ifa, std::ofstream &ofs_running)
         }
     }
 
-    if (PARAM.globalv.rpa_setorb)
-    {
-        if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "ABFS_ORBITAL"))
-        {
-            std::cout << "RPA_EXX_LCAO read abfs_orb!!!" << std::endl;
-            GlobalV::rpa_orbitals.resize(ntype);
-            for (int i = 0; i < ntype; i++)
-            {
-                ifa >> GlobalV::rpa_orbitals[i];
-            }
-        }
-    }
 #endif // __EXX
 #endif // __MPI
 #endif // __LCAO
@@ -640,11 +628,9 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
 
                     if(PARAM.inp.nspin==4)
                     {
-                        if(GlobalV::NONCOLIN)
+                        if(PARAM.inp.noncolin)
                         {
                             //if magnetization only along z-axis, default settings are DOMAG_Z=true and DOMAG=false
-                            GlobalV::DOMAG_Z = false;
-                            GlobalV::DOMAG = true;
                             if(input_angle_mag)
                             {
                                 atoms[it].m_loc_[ia].z = atoms[it].mag[ia] *