Refactor : merging the refactor in lcao and pw

.github/workflows/test.yml

@@ -8,7 +8,7 @@ on:
       - ABACUS_2.2.0_beta
       - deepks
       - planewave
-      - TDDFT
+      - pw_refactor
 jobs:
 
   test:

CMakeLists.txt

@@ -272,6 +272,7 @@ target_link_libraries(${ABACUS_BIN_NAME}
     pw
     ri
     driver
+    en_solver
     -lm
 )
 

doc/input-main.md

@@ -992,7 +992,7 @@ This part of variables are relevant when using hybrid functionals
 
 - exx_hybrid_type<a id="exx-hybrid-type"></a>
     - *Type*: String
-    - *Description*: Type of hybrid functional used. Options are "hf" (pure Hartree-Fock), "pbe0"(PBE0), "hse" (Note: in order to use HSE functional, LIBXC is required).
+    - *Description*: Type of hybrid functional used. Options are "hf" (pure Hartree-Fock), "pbe0"(PBE0), "hse" (Note: in order to use HSE functional, LIBXC is required). Note also that HSE has been tested while PBE0 has NOT been fully tested yet, and the maxmum parallel cpus for running exx is Nx(N+1)/2, with N being the number of atoms.
 
 
         If set to "no", then no hybrid functional is used (i.e.,Fock exchange is not included.)
@@ -1023,7 +1023,7 @@ adial integration for pseudopotentials, in Bohr.
 
 - exx_pca_threshold<a id="exx-pca-threshold"></a>
     - *Type*: Real
-    - *Description*: To accelerate the evaluation of four-center integrals (ik|jl), the product of atomic orbitals are expanded in the basis of auxiliary basis functions (ABF): &phi;<sub>i</sub>&phi;<sub>j</sub>~C<sup>k</sup><sub>ij</sub>P<sub>k</sub>. The size of the ABF (i.e. number of P<sub>k</sub>) is reduced using principal component analysis. When a large PCA threshold is used, the number of ABF will be reduced, hence the calculations becomes faster. However this comes at the cost of computational accuracy. A relatively safe choice of the value is 1d-3.
+    - *Description*: To accelerate the evaluation of four-center integrals (ik|jl), the product of atomic orbitals are expanded in the basis of auxiliary basis functions (ABF): &phi;<sub>i</sub>&phi;<sub>j</sub>~C<sup>k</sup><sub>ij</sub>P<sub>k</sub>. The size of the ABF (i.e. number of P<sub>k</sub>) is reduced using principal component analysis. When a large PCA threshold is used, the number of ABF will be reduced, hence the calculations becomes faster. However this comes at the cost of computational accuracy. A relatively safe choice of the value is 1d-4.
     - *Default*: 0
 
     [back to top](#input-file)
@@ -1044,21 +1044,21 @@ adial integration for pseudopotentials, in Bohr.
 
 - exx_dm_threshold<a id="exx-dm-threshold"></a>
     - *Type*: Real
-    - *Description*: The Fock exchange can be expressed as &Sigma;<sub>k,l</sub>(ik|jl)D<sub>kl</sub> where D is the density matrix. Smaller values of the density matrix can be truncated to accelerate calculation. The larger the threshold is, the faster the calculation and the lower the accuracy. A relatively safe choice of the value is 1d-3.
+    - *Description*: The Fock exchange can be expressed as &Sigma;<sub>k,l</sub>(ik|jl)D<sub>kl</sub> where D is the density matrix. Smaller values of the density matrix can be truncated to accelerate calculation. The larger the threshold is, the faster the calculation and the lower the accuracy. A relatively safe choice of the value is 1d-4.
     - *Default*: 0
 
     [back to top](#input-file)
 
 - exx_schwarz_threshold<a id="exx-schwarz-threshold"></a>
     - *Type*: Real
-    - *Description*: In practice the four-center integrals are sparse, and using Cauchy-Schwartz inequality, we can find an upper bound of each integral before carrying out explicit evaluations. Those that are smaller than exx_schwarz_threshold will be truncated. The larger the threshold is, the faster the calculation and the lower the accuracy. A relatively safe choice of the value is 1d-4.
+    - *Description*: In practice the four-center integrals are sparse, and using Cauchy-Schwartz inequality, we can find an upper bound of each integral before carrying out explicit evaluations. Those that are smaller than exx_schwarz_threshold will be truncated. The larger the threshold is, the faster the calculation and the lower the accuracy. A relatively safe choice of the value is 1d-5.
     - *Default*: 0
 
     [back to top](#input-file)
 
 - exx_cauchy_threshold<a id="exx-cauchy-threshold"></a>
     - *Type*: Real
-    - *Description*: In practice the Fock exchange matrix is sparse, and using Cauchy-Schwartz inequality, we can find an upper bound of each matrix element before carrying out explicit evaluations. Those that are smaller than exx_cauchy_threshold will be truncated. The larger the threshold is, the faster the calculation and the lower the accuracy. A relatively safe choice of the value is 1d-6.
+    - *Description*: In practice the Fock exchange matrix is sparse, and using Cauchy-Schwartz inequality, we can find an upper bound of each matrix element before carrying out explicit evaluations. Those that are smaller than exx_cauchy_threshold will be truncated. The larger the threshold is, the faster the calculation and the lower the accuracy. A relatively safe choice of the value is 1d-7.
     - *Default*: 0
 
     [back to top](#input-file)

source/CMakeLists.txt

@@ -5,6 +5,7 @@ add_subdirectory(module_neighbor)
 add_subdirectory(module_orbital)
 add_subdirectory(module_md)
 add_subdirectory(module_deepks)
+add_subdirectory(module_ensolver)
 add_subdirectory(src_io)
 add_subdirectory(src_ions)
 add_subdirectory(src_lcao)

source/Makefile

@@ -25,6 +25,8 @@ VPATH=./src_global\
 :./src_ri\
 :./\
 
+include module_ensolver/Makefile.ensolver
+
 #==========================
 # Define HONG
 #==========================

source/Makefile.Objects

@@ -176,7 +176,6 @@ FORCE_k.o\
 parallel_orbitals.o \
 global_fp.o \
 pdiag_double.o \
-pdiag_basic.o \
 pdiag_common.o \
 diag_scalapack_gvx.o \
 subgrid_oper.o \

source/driver.cpp

@@ -85,22 +85,31 @@ void Driver::reading(void)
 void Driver::atomic_world(void)
 {
 	ModuleBase::TITLE("Driver","atomic_world");
-
 	//--------------------------------------------------
 	// choose basis sets:
 	// pw: plane wave basis set
 	// lcao_in_pw: LCAO expaned by plane wave basis set
 	// lcao: linear combination of atomic orbitals
 	//--------------------------------------------------
+	string use_ensol;
+	ModuleEnSover::En_Solver *p_ensolver;
 	if(GlobalV::BASIS_TYPE=="pw" || GlobalV::BASIS_TYPE=="lcao_in_pw")
 	{
-		Run_pw::plane_wave_line();
+		use_ensol = "ksdft_pw"; 
+		//We set it temporarily
+		//Finally, we have ksdft_pw, ksdft_lcao, sdft_pw, ofdft, lj, eam, etc.
+		ModuleEnSover::init_esolver(p_ensolver, use_ensol);
+		Run_pw::plane_wave_line(p_ensolver);
+		ModuleEnSover::clean_esolver(p_ensolver);
 	}
 #ifdef __LCAO
 	else if(GlobalV::BASIS_TYPE=="lcao")
-	{
-		Run_lcao::lcao_line();
-	}
+    {
+        use_ensol = "ksdft_lcao";
+        ModuleEnSover::init_esolver(p_ensolver, use_ensol);
+        Run_lcao::lcao_line(p_ensolver);
+        ModuleEnSover::clean_esolver(p_ensolver);
+    }
 #endif
 
 	ModuleBase::timer::finish( GlobalV::ofs_running );

source/module_base/intarray.h

@@ -5,74 +5,152 @@
 #ifndef INTARRAY_H
 #define INTARRAY_H
 
-#include <iostream>
+#include <cassert>
 #include <fstream>
 #include <iomanip>
-#include <cassert>
+#include <iostream>
 
 #ifdef _MCD_CHECK
 //#include "./src_parallel/mcd.h"
 #endif
 
 namespace ModuleBase
 {
-
+/**
+ * @brief Integer array
+ *
+ */
 class IntArray
 {
-public:
-	int *ptr;
-
-	// Constructors for different dimesnions
-	IntArray(const int d1 = 1, const int d2 = 1);
-	IntArray(const int d1, const int d2,const int d3);
-	IntArray(const int d1, const int d2,const int d3,const int d4);
-	IntArray(const int d1, const int d2,const int d3,const int d4,const int d5);
-	IntArray(const int d1, const int d2,const int d3,const int d4,const int d5,const int d6);
-
-	~IntArray();
-
-	void create(const int d1, const int d2);
-	void create(const int d1, const int d2, const int d3);
-	void create(const int d1, const int d2, const int d3, const int d4);
-	void create(const int d1, const int d2, const int d3, const int d4, const int d5);
-	void create(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
-
-	const IntArray &operator=(const IntArray &right);
-	const IntArray &operator=(const int &right);
-
-	int &operator()(const int d1, const int d2);
-	int &operator()(const int d1, const int d2, const int d3);
-	int &operator()(const int d1, const int d2, const int d3,const int d4);
-	int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5);
-	int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
-
-	const int &operator()(const int d1,const int d2)const;
-	const int &operator()(const int d1,const int d2,const int d3)const;
-	const int &operator()(const int d1,const int d2,const int d3,const int d4)const;
-	const int &operator()(const int d1,const int d2,const int d3,const int d4, const int d5)const;
-	const int &operator()(const int d1,const int d2,const int d3,const int d4, const int d5, const int d6)const;
-
-	void zero_out(void);
-
-	int getSize() const{ return size;}
-	int getDim() const{ return dim;}
-	int getBound1() const{ return bound1;}
-	int getBound2() const{ return bound2;}
-	int getBound3() const{ return bound3;}
-	int getBound4() const { return bound4;}
-	int getBound5() const { return bound5;}
-	int getBound6() const { return bound6;}
-
-	static int getArrayCount(void)
-	{ return arrayCount;}
-
-private:
-	int size;
-	int dim;
-	int bound1, bound2, bound3, bound4, bound5, bound6;
-	static int arrayCount;
-	void freemem();
+  public:
+    int *ptr;
+
+    /**
+     * @brief Construct a new Int Array object
+     *
+     * @param d1 The first dimension size
+     * @param d2 The second dimension size
+     */
+    IntArray(const int d1 = 1, const int d2 = 1);
+    IntArray(const int d1, const int d2, const int d3);
+    IntArray(const int d1, const int d2, const int d3, const int d4);
+    IntArray(const int d1, const int d2, const int d3, const int d4, const int d5);
+    IntArray(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
+
+    ~IntArray();
+
+    /**
+     * @brief Create integer arrays
+     *
+     * @param[in] d1
+     * @param[in] d2
+     */
+    void create(const int d1, const int d2);
+    void create(const int d1, const int d2, const int d3);
+    void create(const int d1, const int d2, const int d3, const int d4);
+    void create(const int d1, const int d2, const int d3, const int d4, const int d5);
+    void create(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
+
+    /**
+     * @brief Equal an IntArray to another one
+     *
+     * @param right
+     * @return const IntArray&
+     */
+    const IntArray &operator=(const IntArray &right);
+
+    /**
+     * @brief Equal all elements of an IntArray to an
+     * integer
+     *
+     * @param right
+     * @return const IntArray&
+     */
+    const IntArray &operator=(const int &right);
+
+    /**
+     * @brief Access elements by using operator "()"
+     *
+     * @param d1
+     * @param d2
+     * @return int&
+     */
+    int &operator()(const int d1, const int d2);
+    int &operator()(const int d1, const int d2, const int d3);
+    int &operator()(const int d1, const int d2, const int d3, const int d4);
+    int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5);
+    int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6);
+
+    /**
+     * @brief Access elements by using "()" through pointer
+     * without changing its elements
+     *
+     * @param d1
+     * @param d2
+     * @return const int&
+     */
+    const int &operator()(const int d1, const int d2) const;
+    const int &operator()(const int d1, const int d2, const int d3) const;
+    const int &operator()(const int d1, const int d2, const int d3, const int d4) const;
+    const int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5) const;
+    const int &operator()(const int d1, const int d2, const int d3, const int d4, const int d5, const int d6) const;
+
+    /**
+     * @brief Set all elements of an IntArray to zero
+     *
+     */
+    void zero_out(void);
+
+    int getSize() const
+    {
+        return size;
+    }
+    int getDim() const
+    {
+        return dim;
+    }
+    int getBound1() const
+    {
+        return bound1;
+    }
+    int getBound2() const
+    {
+        return bound2;
+    }
+    int getBound3() const
+    {
+        return bound3;
+    }
+    int getBound4() const
+    {
+        return bound4;
+    }
+    int getBound5() const
+    {
+        return bound5;
+    }
+    int getBound6() const
+    {
+        return bound6;
+    }
+
+    /**
+     * @brief Get the Array Count object
+     *
+     * @return int
+     */
+    static int getArrayCount(void)
+    {
+        return arrayCount;
+    }
+
+  private:
+    int size;
+    int dim;
+    int bound1, bound2, bound3, bound4, bound5, bound6;
+    static int arrayCount;
+    void freemem();
 };
-}
+} // namespace ModuleBase
 
-#endif	// IntArray class
+#endif // IntArray class

source/module_base/math_bspline.cpp

@@ -32,18 +32,12 @@ namespace ModuleBase
         }
     }
 
-    void Bspline::cleanp()
-    {
-        delete[] bezier;
-        bezier = NULL;
-    }
-
     double Bspline::bezier_ele(int n)
     {
         return this->bezier[n];
     }
 
-    void Bspline::getbslpine(double x)
+    void Bspline::getbspline(double x)
     {
         bezier[0] = 1.0;
         for(int k = 1 ; k <= norder ; ++k)