Skip to content

Commit 716e3fd

Browse files
hongriTianqihongriTianqi
authored
Merge branch 'develop' into unittest
2 parents 96a721b + be97d86 commit 716e3fd

File tree

6 files changed

+429
-21
lines changed

6 files changed

+429
-21
lines changed

‎doc/input-main.md‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -992,7 +992,7 @@ This part of variables are relevant when using hybrid functionals
992992
993993
- exx_hybrid_type<a id="exx-hybrid-type"></a>
994994
- *Type*: String
995-
- *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).
995+
- *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.
996996
997997
998998
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.
10231023
10241024
- exx_pca_threshold<a id="exx-pca-threshold"></a>
10251025
- *Type*: Real
1026-
- *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.
1026+
- *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.
10271027
- *Default*: 0
10281028
10291029
[back to top](#input-file)
@@ -1044,21 +1044,21 @@ adial integration for pseudopotentials, in Bohr.
10441044
10451045
- exx_dm_threshold<a id="exx-dm-threshold"></a>
10461046
- *Type*: Real
1047-
- *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.
1047+
- *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.
10481048
- *Default*: 0
10491049
10501050
[back to top](#input-file)
10511051
10521052
- exx_schwarz_threshold<a id="exx-schwarz-threshold"></a>
10531053
- *Type*: Real
1054-
- *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.
1054+
- *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.
10551055
- *Default*: 0
10561056
10571057
[back to top](#input-file)
10581058
10591059
- exx_cauchy_threshold<a id="exx-cauchy-threshold"></a>
10601060
- *Type*: Real
1061-
- *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.
1061+
- *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.
10621062
- *Default*: 0
10631063
10641064
[back to top](#input-file)

‎source/module_base/math_polyint.h‎

Lines changed: 61 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,19 @@ class PolyInt
1818
//========================================================
1919
// Polynomial_Interpolation
2020
//========================================================
21+
22+
/**
23+
* @brief Lagrange interpolation
24+
*
25+
* @param table [in] three dimension matrix, the data in 3rd dimension is used to do prediction
26+
* @param dim1 [in] index of 1st dimension of table/y
27+
* @param dim2 [in] index of 2nd dimension of table/y
28+
* @param y [out] three dimension matrix to store the predicted value
29+
* @param dim_y [in] index of 3rd dimension of y to store predicted value
30+
* @param table_length [in] length of 3rd dimension of table
31+
* @param table_interval [in] interval of 3rd dimension of table
32+
* @param x [in] the position in 3rd dimension to be predicted
33+
*/
2134
static void Polynomial_Interpolation
2235
(
2336
const ModuleBase::realArray &table,
@@ -30,41 +43,84 @@ class PolyInt
3043
const double &x
3144
);
3245

46+
/**
47+
* @brief Lagrange interpolation
48+
*
49+
* @param table [in] three dimension matrix, the data in 3rd dimension is used to do prediction
50+
* @param dim1 [in] index of 1st dimension of table
51+
* @param dim2 [in] index of 2nd dimension of table
52+
* @param table_length [in] length of 3rd dimension of table
53+
* @param table_interval [in] interval of 3rd dimension of table
54+
* @param x [in] the position in 3rd dimension to be predicted
55+
* @return double the predicted value
56+
*/
3357
static double Polynomial_Interpolation
3458
(
3559
const ModuleBase::realArray &table,
3660
const int &dim1,
3761
const int &dim2,
3862
const int &table_length,
3963
const double &table_interval,
40-
const double &x // input value
64+
const double &x
4165
);
4266

43-
static double Polynomial_Interpolation // pengfei Li 2018-3-23
67+
/**
68+
* @brief Lagrange interpolation
69+
*
70+
* @param table [in] four dimension matrix, the data in 4th dimension is used to do prediction
71+
* @param dim1 [in] index of 1st dimension of table
72+
* @param dim2 [in] index of 2nd dimension of table
73+
* @param dim3 [in] index of 3rd dimension of table
74+
* @param table_length [in] length of 4th dimension of table
75+
* @param table_interval [in] interval of 4th dimension of table
76+
* @param x [in] the position in 4th dimension to be predicted
77+
* @return double the predicted value
78+
* @author pengfei Li
79+
* @date 2018-3-23
80+
*/
81+
static double Polynomial_Interpolation
4482
(
4583
const ModuleBase::realArray &table,
4684
const int &dim1,
4785
const int &dim2,
4886
const int &dim3,
4987
const int &table_length,
5088
const double &table_interval,
51-
const double &x // input value
89+
const double &x
5290
);
5391

92+
/**
93+
* @brief Lagrange interpolation
94+
*
95+
* @param table [in] the data used to do prediction
96+
* @param table_length [in] length of table
97+
* @param table_interval [in] interval of table
98+
* @param x [in] the position to be predicted
99+
* @return double the predicted value
100+
*/
54101
static double Polynomial_Interpolation
55102
(
56103
const double *table,
57104
const int &table_length,
58105
const double &table_interval,
59-
const double &x // input value
106+
const double &x
60107
);
61108

109+
/**
110+
* @brief Lagrange interpolation
111+
*
112+
* @param xpoint [in] array of postion
113+
* @param ypoint [in] array of data to do prediction
114+
* @param table_length [in] length of xpoint
115+
* @param x [in] position to be predicted
116+
* @return double predicted value
117+
*/
62118
static double Polynomial_Interpolation_xy
63119
(
64120
const double *xpoint,
65121
const double *ypoint,
66122
const int table_length,
67-
const double &x // input value
123+
const double &x
68124
);
69125

70126
};

‎source/module_base/math_sphbes.h‎

Lines changed: 38 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -14,26 +14,53 @@ class Sphbes
1414
Sphbes();
1515
~Sphbes();
1616

17+
/**
18+
* @brief spherical bessel
19+
*
20+
* @param msh [in] number of grid points
21+
* @param r [in] radial grid
22+
* @param q [in] k_radial
23+
* @param l [in] angular momentum
24+
* @param jl [out] jl spherical bessel function
25+
*/
1726
static void Spherical_Bessel
1827
(
19-
const int &msh, //number of grid points
20-
const double *r,//radial grid
21-
const double &q, //
22-
const int &l, //angular momentum
23-
double *jl //jl(1:msh) = j_l(q*r(i)),spherical bessel function
28+
const int &msh,
29+
const double *r,
30+
const double &q,
31+
const int &l,
32+
double *jl
2433
);
2534

35+
/**
36+
* @brief spherical bessel
37+
*
38+
* @param msh [in] number of grid points
39+
* @param r [in] radial grid
40+
* @param q [in] k_radial
41+
* @param l [in] angular momentum
42+
* @param jl [out] jl spherical bessel function
43+
* @param sjp [out] sjp[i] is assigned to be 1.0. i < msh.
44+
*/
2645
static void Spherical_Bessel
2746
(
28-
const int &msh, //number of grid points
29-
const double *r,//radial grid
30-
const double &q, //
31-
const int &l, //angular momentum
32-
double *sj, //jl(1:msh) = j_l(q*r(i)),spherical bessel function
47+
const int &msh,
48+
const double *r,
49+
const double &q,
50+
const int &l,
51+
double *sj,
3352
double *sjp
3453
);
3554

36-
55+
/**
56+
* @brief return num eigenvalues of spherical bessel function
57+
*
58+
* @param num [in] the number of eigenvalues
59+
* @param l [in] angular number
60+
* @param epsilon [in] the accuracy
61+
* @param eigenvalue [out] the calculated eigenvalues
62+
* @param rcut [in] the cutoff the radial function
63+
*/
3764
static void Spherical_Bessel_Roots
3865
(
3966
const int &num,

‎source/module_base/test/CMakeLists.txt‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,10 @@ AddTest(
4747
TARGET base_sph_bessel_recursive
4848
SOURCES sph_bessel_recursive_test.cpp ../sph_bessel_recursive-d1.cpp ../sph_bessel_recursive-d2.cpp
4949
)
50+
AddTest(
51+
TARGET base_math_sphbes
52+
SOURCES math_sphbes_test.cpp ../math_sphbes.cpp ../timer.cpp
53+
)
5054
AddTest(
5155
TARGET base_realarray
5256
SOURCES realarray_test.cpp ../realarray.cpp
@@ -64,6 +68,10 @@ AddTest(
6468
LIBS ${math_libs}
6569
SOURCES mathzone_test.cpp ../mathzone.cpp ../matrix3.cpp ../matrix.cpp ../tool_quit.cpp ../global_variable.cpp ../global_file.cpp ../memory.cpp ../timer.cpp
6670
)
71+
AddTest(
72+
TARGET base_math_polyint
73+
SOURCES math_polyint_test.cpp ../math_polyint.cpp ../realarray.cpp ../timer.cpp
74+
)
6775
AddTest(
6876
TARGET base_mathzone_add1
6977
LIBS ${math_libs}

‎source/module_base/test/math_polyint_test.cpp‎

Lines changed: 136 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,136 @@
1+
#include"../math_polyint.h"
2+
#include"gtest/gtest.h"
3+
#include"../realarray.h"
4+
#include<math.h>
5+
6+
#define doublethreshold 1e-9
7+
8+
/************************************************
9+
* unit test of class PolyInt
10+
***********************************************/
11+
12+
/**
13+
* This unit test is to verify the accuracy of
14+
* interpolation method on the function sin(x)/x
15+
* with a interval of 0.01.
16+
* sin(x)/x is one of the solution of spherical bessel
17+
* function when l=0.
18+
*
19+
* - Tested function:
20+
* - 4 types of Polynomial_Interpolation
21+
* - Polynomial_Interpolation_xy
22+
*/
23+
24+
25+
class bessell0 : public testing::Test
26+
{
27+
protected:
28+
29+
int TableLength = 400;
30+
double interval = 0.01;
31+
ModuleBase::realArray table3,table4;
32+
ModuleBase::realArray y3;
33+
double *tablex;
34+
double *tabley;
35+
36+
double sinc(double x) {return sin(x)/x;}
37+
38+
void SetUp()
39+
{
40+
tablex = new double[TableLength];
41+
tabley = new double[TableLength];
42+
table3.create(1,1,TableLength);
43+
table4.create(1,1,1,TableLength);
44+
y3.create(1,1,TableLength);
45+
46+
for(int i=1;i<TableLength;++i)
47+
{
48+
table3(0,0,i) = sinc(i * interval);
49+
table4(0,0,0,i) = sinc(i * interval);
50+
tablex[i] = i * interval;
51+
tabley[i] = sinc(i * interval);
52+
}
53+
}
54+
55+
void TearDown()
56+
{
57+
delete [] tablex;
58+
delete [] tabley;
59+
}
60+
};
61+
62+
63+
TEST_F(bessell0,PolynomialInterpolationThreeDimensionY)
64+
{
65+
ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,1,TableLength,interval,0.1);
66+
ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,2,TableLength,interval,1.005);
67+
ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,3,TableLength,interval,2.005);
68+
ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,4,TableLength,interval,3.005);
69+
ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,5,TableLength,interval,3.505);
70+
71+
EXPECT_NEAR(y3(0,0,1),sinc(0.1),doublethreshold);
72+
EXPECT_NEAR(y3(0,0,2),sinc(1.005),doublethreshold);
73+
EXPECT_NEAR(y3(0,0,3),sinc(2.005),doublethreshold);
74+
EXPECT_NEAR(y3(0,0,4),sinc(3.005),doublethreshold);
75+
EXPECT_NEAR(y3(0,0,5),sinc(3.505),doublethreshold);
76+
}
77+
78+
TEST_F(bessell0,PolynomialInterpolationThreeDimension)
79+
{
80+
double y1 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,0.1);
81+
double y2 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,1.005);
82+
double y3 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,2.005);
83+
double y4 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,3.005);
84+
double y5 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,3.505);
85+
86+
EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
87+
EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
88+
EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
89+
EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
90+
EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
91+
}
92+
93+
TEST_F(bessell0,PolynomialInterpolationFourDimension)
94+
{
95+
double y1 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,0.1);
96+
double y2 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,1.005);
97+
double y3 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,2.005);
98+
double y4 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,3.005);
99+
double y5 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,3.505);
100+
101+
EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
102+
EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
103+
EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
104+
EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
105+
EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
106+
}
107+
108+
TEST_F(bessell0,PolynomialInterpolation)
109+
{
110+
double y1 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,0.1);
111+
double y2 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,1.005);
112+
double y3 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,2.005);
113+
double y4 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,3.005);
114+
double y5 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,3.505);
115+
116+
EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
117+
EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
118+
EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
119+
EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
120+
EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
121+
}
122+
123+
TEST_F(bessell0,PolynomialInterpolationXY)
124+
{
125+
double y1 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,0.1);
126+
double y2 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,1.005);
127+
double y3 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,2.005);
128+
double y4 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,3.005);
129+
double y5 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,3.505);
130+
131+
EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
132+
EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
133+
EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
134+
EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
135+
EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
136+
}

0 commit comments

Comments
 (0)