diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 41318c963ec..e7cdd42dc31 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -8,7 +8,7 @@ on:
- ABACUS_2.2.0_beta
- deepks
- planewave
- - TDDFT
+ - pw_refactor
jobs:
test:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1d4415a82d4..8716afbc8c7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -272,6 +272,7 @@ target_link_libraries(${ABACUS_BIN_NAME}
pw
ri
driver
+ en_solver
-lm
)
diff --git a/doc/input-main.md b/doc/input-main.md
index 654ce8714de..8e87f6175b5 100644
--- a/doc/input-main.md
+++ b/doc/input-main.md
@@ -992,7 +992,7 @@ This part of variables are relevant when using hybrid functionals
- exx_hybrid_type
- *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
- *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): φiφj~CkijPk. The size of the ABF (i.e. number of Pk) 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): φiφj~CkijPk. The size of the ABF (i.e. number of Pk) 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
- *Type*: Real
- - *Description*: The Fock exchange can be expressed as Σk,l(ik|jl)Dkl 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 Σk,l(ik|jl)Dkl 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
- *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
- *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)
diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt
index e65955d5989..fb72efffb62 100644
--- a/source/CMakeLists.txt
+++ b/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)
diff --git a/source/Makefile b/source/Makefile
index b3bcc1d2218..eab82838029 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -25,6 +25,8 @@ VPATH=./src_global\
:./src_ri\
:./\
+include module_ensolver/Makefile.ensolver
+
#==========================
# Define HONG
#==========================
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index 44a72c0d0e2..f486266f201 100644
--- a/source/Makefile.Objects
+++ b/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 \
diff --git a/source/driver.cpp b/source/driver.cpp
index 6128759c1d9..42412311a2d 100644
--- a/source/driver.cpp
+++ b/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 );
diff --git a/source/module_base/intarray.h b/source/module_base/intarray.h
index f8db1250753..64cacf014b7 100644
--- a/source/module_base/intarray.h
+++ b/source/module_base/intarray.h
@@ -5,10 +5,10 @@
#ifndef INTARRAY_H
#define INTARRAY_H
-#include
+#include
#include
#include
-#include
+#include
#ifdef _MCD_CHECK
//#include "./src_parallel/mcd.h"
@@ -16,63 +16,141 @@
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
diff --git a/source/module_base/math_bspline.cpp b/source/module_base/math_bspline.cpp
index 378fcd1ed5e..63d95e3f97b 100644
--- a/source/module_base/math_bspline.cpp
+++ b/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)
diff --git a/source/module_base/math_bspline.h b/source/module_base/math_bspline.h
index 2eb72d8880e..86155bb79dc 100644
--- a/source/module_base/math_bspline.h
+++ b/source/module_base/math_bspline.h
@@ -4,57 +4,54 @@
namespace ModuleBase
{
-
-//
-//DESCRIPTION:
-// A class to treat Cardinal B-spline interpolation.
-// qianrui created 2021-09-14
-//MATH:
-// Only uniform nodes are considered: xm-x[m-1]=Dx(>= 0) for control node: X={x0,x1,...,xm};
-// Any function p(x) can be written by
-// p(x)=\sum_i{ci*M_ik(x)} (k->infinity),
-// where M_ik is the i-th k-order Cardinal B-spline base function
-// and ci is undetermined coefficient.
-// M_i0 = H(x-xi)-H(x-x[i+1]), H(x): step function
-// x-xi x[i+k+1]-x
-// M_ik(x)= ---------*M_i(k-1)(x)+ ----------------*M_[i+1][k-1](x) ( xi <= x <= x[i+1] )
-// x[i+k]-xi x[i+k+1]-x[i+1]
-// For uniform nodes: M_[i+1]k(x+Dx)=M_ik(x)
-// If we define Bk[n] stores M_ik(x+n*Dx) for x in (xi,xi+Dx):
-// x+n*Dx-xi xi+(k-n+1)*Dx-x
-// Bk[n] = -----------*B(k-1)[n] + -----------------*B(k-1)[n-1]
-// k*Dx k*Dx
-//USAGE:
-// ModuleBase::Bspline bp;
-// bp.init(10,0.7,2); //Dx = 0.7, xi = 2
-// bp.getbslpine(0.5); //x = 0.5
-// cout<= 0) for control node: X={x0,x1,...,xm};
+ * Any function p(x) can be written by
+ * p(x)=\sum_i{ci*M_ik(x)} (k->infinity),
+ * where M_ik is the i-th k-order Cardinal B-spline base function
+ * and ci is undetermined coefficient.
+ * M_i0 = H(x-xi)-H(x-x[i+1]), H(x): step function
+ * x-xi x[i+k+1]-x
+ * M_ik(x)= ---------*M_i(k-1)(x)+ ----------------*M_[i+1][k-1](x) ( xi <= x <= x[i+1] )
+ * x[i+k]-xi x[i+k+1]-x[i+1]
+ * For uniform nodes: M_[i+1]k(x+Dx)=M_ik(x)
+ * If we define Bk[n] stores M_ik(x+n*Dx) for x in (xi,xi+Dx):
+ * x+n*Dx-xi xi+(k-n+1)*Dx-x
+ * Bk[n] = -----------*B(k-1)[n] + -----------------*B(k-1)[n-1]
+ * k*Dx k*Dx
+ * USAGE:
+ * ModuleBase::Bspline bp;
+ * bp.init(10,0.7,2); //Dx = 0.7, xi = 2
+ * bp.getbslpine(0.5); //x = 0.5
+ * cout<= 0
- double Dx; //Dx: the interval of control node
- double xi; // xi: the starting point
- double *bezier; //bezier[n] = Bk[n]
+ private:
+ int norder; // the order of bezier base; norder >= 0
+ double Dx; // Dx: the interval of control node
+ double xi; // xi: the starting point
+ double *bezier; // bezier[n] = Bk[n]
- public:
- Bspline();
- ~Bspline();
+ public:
+ Bspline();
+ ~Bspline();
- //Init norder, Dx, xi
- void init(int norderin, double Dxin, double xiin);
+ void init(int norderin, double Dxin, double xiin);
- //delete[] bezier
- void cleanp();
+ // Get the result of i-th bezier base functions for different input x+xi+n*Dx.
+ // x should be in [0,Dx]
+ // n-th result is stored in bezier[n];
+ void getbspline(double x);
- //Get the result of i-th bezier base functions for different input x+xi+n*Dx.
- //x should be in [0,Dx]
- //n-th result is stored in bezier[n];
- void getbslpine(double x);
-
- //get the element of bezier
- double bezier_ele(int n);
+ // get the element of bezier
+ double bezier_ele(int n);
};
-}
+} // namespace ModuleBase
#endif
diff --git a/source/module_base/math_polyint.h b/source/module_base/math_polyint.h
index 0cc13f8f831..b15f50678c3 100644
--- a/source/module_base/math_polyint.h
+++ b/source/module_base/math_polyint.h
@@ -18,6 +18,19 @@ class PolyInt
//========================================================
// Polynomial_Interpolation
//========================================================
+
+ /**
+ * @brief Lagrange interpolation
+ *
+ * @param table [in] three dimension matrix, the data in 3rd dimension is used to do prediction
+ * @param dim1 [in] index of 1st dimension of table/y
+ * @param dim2 [in] index of 2nd dimension of table/y
+ * @param y [out] three dimension matrix to store the predicted value
+ * @param dim_y [in] index of 3rd dimension of y to store predicted value
+ * @param table_length [in] length of 3rd dimension of table
+ * @param table_interval [in] interval of 3rd dimension of table
+ * @param x [in] the position in 3rd dimension to be predicted
+ */
static void Polynomial_Interpolation
(
const ModuleBase::realArray &table,
@@ -30,6 +43,17 @@ class PolyInt
const double &x
);
+ /**
+ * @brief Lagrange interpolation
+ *
+ * @param table [in] three dimension matrix, the data in 3rd dimension is used to do prediction
+ * @param dim1 [in] index of 1st dimension of table
+ * @param dim2 [in] index of 2nd dimension of table
+ * @param table_length [in] length of 3rd dimension of table
+ * @param table_interval [in] interval of 3rd dimension of table
+ * @param x [in] the position in 3rd dimension to be predicted
+ * @return double the predicted value
+ */
static double Polynomial_Interpolation
(
const ModuleBase::realArray &table,
@@ -37,10 +61,24 @@ class PolyInt
const int &dim2,
const int &table_length,
const double &table_interval,
- const double &x // input value
+ const double &x
);
- static double Polynomial_Interpolation // pengfei Li 2018-3-23
+ /**
+ * @brief Lagrange interpolation
+ *
+ * @param table [in] four dimension matrix, the data in 4th dimension is used to do prediction
+ * @param dim1 [in] index of 1st dimension of table
+ * @param dim2 [in] index of 2nd dimension of table
+ * @param dim3 [in] index of 3rd dimension of table
+ * @param table_length [in] length of 4th dimension of table
+ * @param table_interval [in] interval of 4th dimension of table
+ * @param x [in] the position in 4th dimension to be predicted
+ * @return double the predicted value
+ * @author pengfei Li
+ * @date 2018-3-23
+ */
+ static double Polynomial_Interpolation
(
const ModuleBase::realArray &table,
const int &dim1,
@@ -48,23 +86,41 @@ class PolyInt
const int &dim3,
const int &table_length,
const double &table_interval,
- const double &x // input value
+ const double &x
);
+ /**
+ * @brief Lagrange interpolation
+ *
+ * @param table [in] the data used to do prediction
+ * @param table_length [in] length of table
+ * @param table_interval [in] interval of table
+ * @param x [in] the position to be predicted
+ * @return double the predicted value
+ */
static double Polynomial_Interpolation
(
const double *table,
const int &table_length,
const double &table_interval,
- const double &x // input value
+ const double &x
);
+ /**
+ * @brief Lagrange interpolation
+ *
+ * @param xpoint [in] array of postion
+ * @param ypoint [in] array of data to do prediction
+ * @param table_length [in] length of xpoint
+ * @param x [in] position to be predicted
+ * @return double predicted value
+ */
static double Polynomial_Interpolation_xy
(
const double *xpoint,
const double *ypoint,
const int table_length,
- const double &x // input value
+ const double &x
);
};
diff --git a/source/module_base/math_sphbes.h b/source/module_base/math_sphbes.h
index f152e52e659..f7257b2dc3e 100644
--- a/source/module_base/math_sphbes.h
+++ b/source/module_base/math_sphbes.h
@@ -14,26 +14,53 @@ class Sphbes
Sphbes();
~Sphbes();
+ /**
+ * @brief spherical bessel
+ *
+ * @param msh [in] number of grid points
+ * @param r [in] radial grid
+ * @param q [in] k_radial
+ * @param l [in] angular momentum
+ * @param jl [out] jl spherical bessel function
+ */
static void Spherical_Bessel
(
- const int &msh, //number of grid points
- const double *r,//radial grid
- const double &q, //
- const int &l, //angular momentum
- double *jl //jl(1:msh) = j_l(q*r(i)),spherical bessel function
+ const int &msh,
+ const double *r,
+ const double &q,
+ const int &l,
+ double *jl
);
+ /**
+ * @brief spherical bessel
+ *
+ * @param msh [in] number of grid points
+ * @param r [in] radial grid
+ * @param q [in] k_radial
+ * @param l [in] angular momentum
+ * @param jl [out] jl spherical bessel function
+ * @param sjp [out] sjp[i] is assigned to be 1.0. i < msh.
+ */
static void Spherical_Bessel
(
- const int &msh, //number of grid points
- const double *r,//radial grid
- const double &q, //
- const int &l, //angular momentum
- double *sj, //jl(1:msh) = j_l(q*r(i)),spherical bessel function
+ const int &msh,
+ const double *r,
+ const double &q,
+ const int &l,
+ double *sj,
double *sjp
);
-
+ /**
+ * @brief return num eigenvalues of spherical bessel function
+ *
+ * @param num [in] the number of eigenvalues
+ * @param l [in] angular number
+ * @param epsilon [in] the accuracy
+ * @param eigenvalue [out] the calculated eigenvalues
+ * @param rcut [in] the cutoff the radial function
+ */
static void Spherical_Bessel_Roots
(
const int &num,
diff --git a/source/module_base/math_ylmreal.h b/source/module_base/math_ylmreal.h
index 362f2dfdd82..c051ef95ca3 100644
--- a/source/module_base/math_ylmreal.h
+++ b/source/module_base/math_ylmreal.h
@@ -14,29 +14,54 @@ class YlmReal
YlmReal();
~YlmReal();
+ /**
+ * @brief spherical harmonic function (real form) an array of vectors
+ *
+ * @param lmax2 [in] lmax2 = (lmax + 1)^2 ; lmax = angular quantum number
+ * @param ng [in] the number of vectors
+ * @param g [in] an array of vectors
+ * @param ylm [out] Ylm; column index represent vector, row index represent Y00, Y10, Y11, Y1-1, Y20,Y21,Y2-1,Y22.Y2-2,...;
+ */
static void Ylm_Real
(
- const int lmax2, // lmax2 = (lmax+1)^2
- const int ng, //
- const ModuleBase::Vector3 *g, // g_cartesian_vec(x,y,z)
- matrix &ylm // output
+ const int lmax2,
+ const int ng,
+ const ModuleBase::Vector3 *g,
+ matrix &ylm
);
+ /**
+ * @brief spherical harmonic function (Herglotz generating form) of an array of vectors
+ *
+ * @param lmax2 [in] lmax2 = (lmax + 1)^2 ; lmax = angular quantum number
+ * @param ng [in] the number of vectors
+ * @param g [in] an array of vectors
+ * @param ylm [out] Ylm; column index represent vector, row index represent Y00, Y10, Y11, Y1-1, Y20,Y21,Y2-1,Y22.Y2-2,...;
+ */
static void Ylm_Real2
(
- const int lmax2, // lmax2 = (lmax+1)^2
- const int ng, //
- const ModuleBase::Vector3 *g, // g_cartesian_vec(x,y,z)
- matrix &ylm // output
+ const int lmax2,
+ const int ng,
+ const ModuleBase::Vector3 *g,
+ matrix &ylm
);
+ /**
+ * @brief spherical harmonic function (Herglotz generating form) of a vector
+ *
+ * @param lmax [in] maximum angular quantum number
+ * @param x [in] x part of the vector
+ * @param y [in] y part of the vector
+ * @param z [in] z part of the vector
+ * @param rly [in] Ylm, Y00, Y10, Y11, Y1-1, Y20,Y21,Y2-1,Y22.Y2-2,...
+ */
static void rlylm
(
const int lmax,
const double& x,
const double& y,
- const double& z, // g_cartesian_vec(x,y,z)
- double* rly // output
+ const double& z,
+ double* rly
);
private:
diff --git a/source/module_base/mathzone.h b/source/module_base/mathzone.h
index 99622ac29a1..6e5180afd8c 100644
--- a/source/module_base/mathzone.h
+++ b/source/module_base/mathzone.h
@@ -1,69 +1,82 @@
#ifndef MATHZONE_H
#define MATHZONE_H
-#include "realarray.h"
-#include "matrix3.h"
#include "global_function.h"
-#include
-#include