Skip to content

Commit 129e193

Browse files
pxlxingliangpxlxingliang
authored
Merge pull request #727 from pxlxingliang/develop
test:add the unittest and comments of math_ylmreal.h
2 parents be97d86 + a3cd7d8 commit 129e193

File tree

3 files changed

+244
-10
lines changed

3 files changed

+244
-10
lines changed

‎source/module_base/math_ylmreal.h‎

Lines changed: 35 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -14,29 +14,54 @@ class YlmReal
1414
YlmReal();
1515
~YlmReal();
1616

17+
/**
18+
* @brief spherical harmonic function (real form) an array of vectors
19+
*
20+
* @param lmax2 [in] lmax2 = (lmax + 1)^2 ; lmax = angular quantum number
21+
* @param ng [in] the number of vectors
22+
* @param g [in] an array of vectors
23+
* @param ylm [out] Ylm; column index represent vector, row index represent Y00, Y10, Y11, Y1-1, Y20,Y21,Y2-1,Y22.Y2-2,...;
24+
*/
1725
static void Ylm_Real
1826
(
19-
const int lmax2, // lmax2 = (lmax+1)^2
20-
const int ng, //
21-
const ModuleBase::Vector3<double> *g, // g_cartesian_vec(x,y,z)
22-
matrix &ylm // output
27+
const int lmax2,
28+
const int ng,
29+
const ModuleBase::Vector3<double> *g,
30+
matrix &ylm
2331
);
2432

33+
/**
34+
* @brief spherical harmonic function (Herglotz generating form) of an array of vectors
35+
*
36+
* @param lmax2 [in] lmax2 = (lmax + 1)^2 ; lmax = angular quantum number
37+
* @param ng [in] the number of vectors
38+
* @param g [in] an array of vectors
39+
* @param ylm [out] Ylm; column index represent vector, row index represent Y00, Y10, Y11, Y1-1, Y20,Y21,Y2-1,Y22.Y2-2,...;
40+
*/
2541
static void Ylm_Real2
2642
(
27-
const int lmax2, // lmax2 = (lmax+1)^2
28-
const int ng, //
29-
const ModuleBase::Vector3<double> *g, // g_cartesian_vec(x,y,z)
30-
matrix &ylm // output
43+
const int lmax2,
44+
const int ng,
45+
const ModuleBase::Vector3<double> *g,
46+
matrix &ylm
3147
);
3248

49+
/**
50+
* @brief spherical harmonic function (Herglotz generating form) of a vector
51+
*
52+
* @param lmax [in] maximum angular quantum number
53+
* @param x [in] x part of the vector
54+
* @param y [in] y part of the vector
55+
* @param z [in] z part of the vector
56+
* @param rly [in] Ylm, Y00, Y10, Y11, Y1-1, Y20,Y21,Y2-1,Y22.Y2-2,...
57+
*/
3358
static void rlylm
3459
(
3560
const int lmax,
3661
const double& x,
3762
const double& y,
38-
const double& z, // g_cartesian_vec(x,y,z)
39-
double* rly // output
63+
const double& z,
64+
double* rly
4065
);
4166

4267
private:

‎source/module_base/test/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,4 +71,9 @@ AddTest(
7171
AddTest(
7272
TARGET base_math_polyint
7373
SOURCES math_polyint_test.cpp ../math_polyint.cpp ../realarray.cpp ../timer.cpp
74+
)
75+
AddTest(
76+
TARGET base_math_ylmreal
77+
LIBS ${math_libs}
78+
SOURCES math_ylmreal_test.cpp ../math_ylmreal.cpp ../realarray.cpp ../timer.cpp ../matrix.cpp
7479
)

‎source/module_base/test/math_ylmreal_test.cpp‎

Lines changed: 204 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,204 @@
1+
#include"../math_ylmreal.h"
2+
#include"../vector3.h"
3+
#include"../matrix.h"
4+
#include"gtest/gtest.h"
5+
#include<math.h>
6+
7+
#define doublethreshold 1e-12
8+
9+
/************************************************
10+
* unit test of class YlmReal and Ylm
11+
***********************************************/
12+
13+
/**
14+
* For lmax <5 cases, the reference values are calculated by the formula from
15+
* https://formulasearchengine.com/wiki/Table_of_spherical_harmonics. Note, these
16+
* formula lack of the Condon–Shortley phase (-1)^m, and in this unit test, item
17+
* (-1)^m is multiplied.
18+
* For lmax >=5, the reference values are calculated by YlmReal::Ylm_Real.
19+
*
20+
* - Tested functions of class YlmReal
21+
* - Ylm_Real
22+
* - Ylm_Real2
23+
* - rlylm
24+
*/
25+
26+
27+
28+
//mock functions of WARNING_QUIT and WARNING
29+
namespace ModuleBase
30+
{
31+
void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
32+
void WARNING(const std::string &file,const std::string &description) {return ;}
33+
}
34+
35+
36+
class YlmRealTest : public testing::Test
37+
{
38+
protected:
39+
40+
int lmax = 7;
41+
int ng = 4; //test the 4 selected points on the sphere
42+
int nylm ; // total Ylm number;
43+
ModuleBase::matrix ylm;
44+
ModuleBase::Vector3<double> *g;
45+
double *ref;
46+
double *rly;
47+
48+
//Ylm function
49+
//https://formulasearchengine.com/wiki/Table_of_spherical_harmonics
50+
//multipy the Condon–Shortley phase (-1)^m
51+
inline double norm(const double &x, const double &y, const double &z) {return sqrt(x*x + y*y + z*z);}
52+
double y00(const double &x, const double &y, const double &z) {return 1.0/2.0/sqrt(M_PI);}
53+
double y10(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return sqrt(3.0/(4.0*M_PI)) * z / r;}
54+
double y11(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*sqrt(3.0/(4.*M_PI)) * x / r;}
55+
double y1m1(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*sqrt(3./(4.*M_PI)) * y / r;} // y1m1 means Y1,-1
56+
double y20(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 1./4. * sqrt(5./M_PI) * (-1.*x*x - y*y + 2.*z*z) / (r*r);}
57+
double y21(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*1./2. * sqrt(15./M_PI) * (z*x) / (r*r);}
58+
double y2m1(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*1./2. * sqrt(15./M_PI) * (z*y) / (r*r);}
59+
double y22(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 1./4. * sqrt(15./M_PI) * (x*x - y*y) / (r*r);}
60+
double y2m2(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 1./2. * sqrt(15./M_PI) * (x*y) / (r*r);}
61+
double y30(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 1./4. * sqrt(7./M_PI) * z*(2.*z*z-3.*x*x-3.*y*y) / (r*r*r);}
62+
double y31(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*1./4. * sqrt(21./2./M_PI) * x*(4.*z*z-x*x-y*y) / (r*r*r);}
63+
double y3m1(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*1./4. * sqrt(21./2./M_PI) * y*(4.*z*z-x*x-y*y) / (r*r*r);}
64+
double y32(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 1./4. * sqrt(105./M_PI) * (x*x - y*y)*z / (r*r*r);}
65+
double y3m2(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 1./2. * sqrt(105./M_PI) * x*y*z / (r*r*r);}
66+
double y33(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*1./4. * sqrt(35./2./M_PI) * x*(x*x - 3.*y*y) / (r*r*r);}
67+
double y3m3(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*1./4. * sqrt(35./2./M_PI) * y*(3.*x*x - y*y) / (r*r*r);}
68+
double y40(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 3./16.*sqrt(1./M_PI) * (35.*z*z*z*z - 30.*z*z*r*r + 3*r*r*r*r) / (r*r*r*r);}
69+
double y41(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*3./4.*sqrt(5./2./M_PI) * x*z*(7.*z*z - 3*r*r) / (r*r*r*r);}
70+
double y4m1(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*3./4.*sqrt(5./2./M_PI) * y*z*(7.*z*z - 3.*r*r) / (r*r*r*r);}
71+
double y42(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 3./8.*sqrt(5./M_PI) * (x*x-y*y)*(7.*z*z-r*r) / (r*r*r*r);}
72+
double y4m2(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 3./4.*sqrt(5./M_PI) * x*y*(7.*z*z - r*r) / (r*r*r*r);}
73+
double y43(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*3./4.*sqrt(35./2./M_PI) * x*z*(x*x - 3.*y*y) / (r*r*r*r);}
74+
double y4m3(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return -1.0*3./4.*sqrt(35./2./M_PI) * y*z*(3.*x*x - y*y) / (r*r*r*r);}
75+
double y44(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 3./16.*sqrt(35./M_PI) * (x*x*(x*x - 3.*y*y) - y*y*(3.*x*x-y*y)) / (r*r*r*r);}
76+
double y4m4(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 3./4.*sqrt(35./M_PI) * x*y*(x*x - y*y) / (r*r*r*r);}
77+
78+
void SetUp()
79+
{
80+
nylm = (lmax + 1) * (lmax + 1);
81+
ylm.create(nylm,ng);
82+
g = new ModuleBase::Vector3<double>[ng];
83+
g[0].set(1.0,0.0,0.0);
84+
g[1].set(0.0,1.0,0.0);
85+
g[2].set(0.0,0.0,1.0);
86+
g[3].set(-1.0,-1.0,-1.0);
87+
88+
rly = new double[nylm];
89+
ref = new double[nylm*ng]{
90+
y00(g[0].x, g[0].y, g[0].z), y00(g[1].x, g[1].y, g[1].z), y00(g[2].x, g[2].y, g[2].z), y00(g[3].x, g[3].y, g[3].z),
91+
y10(g[0].x, g[0].y, g[0].z), y10(g[1].x, g[1].y, g[1].z), y10(g[2].x, g[2].y, g[2].z), y10(g[3].x, g[3].y, g[3].z),
92+
y11(g[0].x, g[0].y, g[0].z), y11(g[1].x, g[1].y, g[1].z), y11(g[2].x, g[2].y, g[2].z), y11(g[3].x, g[3].y, g[3].z),
93+
y1m1(g[0].x, g[0].y, g[0].z), y1m1(g[1].x, g[1].y, g[1].z), y1m1(g[2].x, g[2].y, g[2].z), y1m1(g[3].x, g[3].y, g[3].z),
94+
y20(g[0].x, g[0].y, g[0].z), y20(g[1].x, g[1].y, g[1].z), y20(g[2].x, g[2].y, g[2].z), y20(g[3].x, g[3].y, g[3].z),
95+
y21(g[0].x, g[0].y, g[0].z), y21(g[1].x, g[1].y, g[1].z), y21(g[2].x, g[2].y, g[2].z), y21(g[3].x, g[3].y, g[3].z),
96+
y2m1(g[0].x, g[0].y, g[0].z), y2m1(g[1].x, g[1].y, g[1].z), y2m1(g[2].x, g[2].y, g[2].z), y2m1(g[3].x, g[3].y, g[3].z),
97+
y22(g[0].x, g[0].y, g[0].z), y22(g[1].x, g[1].y, g[1].z), y22(g[2].x, g[2].y, g[2].z), y22(g[3].x, g[3].y, g[3].z),
98+
y2m2(g[0].x, g[0].y, g[0].z), y2m2(g[1].x, g[1].y, g[1].z), y2m2(g[2].x, g[2].y, g[2].z), y2m2(g[3].x, g[3].y, g[3].z),
99+
y30(g[0].x, g[0].y, g[0].z), y30(g[1].x, g[1].y, g[1].z), y30(g[2].x, g[2].y, g[2].z), y30(g[3].x, g[3].y, g[3].z),
100+
y31(g[0].x, g[0].y, g[0].z), y31(g[1].x, g[1].y, g[1].z), y31(g[2].x, g[2].y, g[2].z), y31(g[3].x, g[3].y, g[3].z),
101+
y3m1(g[0].x, g[0].y, g[0].z), y3m1(g[1].x, g[1].y, g[1].z), y3m1(g[2].x, g[2].y, g[2].z), y3m1(g[3].x, g[3].y, g[3].z),
102+
y32(g[0].x, g[0].y, g[0].z), y32(g[1].x, g[1].y, g[1].z), y32(g[2].x, g[2].y, g[2].z), y32(g[3].x, g[3].y, g[3].z),
103+
y3m2(g[0].x, g[0].y, g[0].z), y3m2(g[1].x, g[1].y, g[1].z), y3m2(g[2].x, g[2].y, g[2].z), y3m2(g[3].x, g[3].y, g[3].z),
104+
y33(g[0].x, g[0].y, g[0].z), y33(g[1].x, g[1].y, g[1].z), y33(g[2].x, g[2].y, g[2].z), y33(g[3].x, g[3].y, g[3].z),
105+
y3m3(g[0].x, g[0].y, g[0].z), y3m3(g[1].x, g[1].y, g[1].z), y3m3(g[2].x, g[2].y, g[2].z), y3m3(g[3].x, g[3].y, g[3].z),
106+
y40(g[0].x, g[0].y, g[0].z), y40(g[1].x, g[1].y, g[1].z), y40(g[2].x, g[2].y, g[2].z), y40(g[3].x, g[3].y, g[3].z),
107+
y41(g[0].x, g[0].y, g[0].z), y41(g[1].x, g[1].y, g[1].z), y41(g[2].x, g[2].y, g[2].z), y41(g[3].x, g[3].y, g[3].z),
108+
y4m1(g[0].x, g[0].y, g[0].z), y4m1(g[1].x, g[1].y, g[1].z), y4m1(g[2].x, g[2].y, g[2].z), y4m1(g[3].x, g[3].y, g[3].z),
109+
y42(g[0].x, g[0].y, g[0].z), y42(g[1].x, g[1].y, g[1].z), y42(g[2].x, g[2].y, g[2].z), y42(g[3].x, g[3].y, g[3].z),
110+
y4m2(g[0].x, g[0].y, g[0].z), y4m2(g[1].x, g[1].y, g[1].z), y4m2(g[2].x, g[2].y, g[2].z), y4m2(g[3].x, g[3].y, g[3].z),
111+
y43(g[0].x, g[0].y, g[0].z), y43(g[1].x, g[1].y, g[1].z), y43(g[2].x, g[2].y, g[2].z), y43(g[3].x, g[3].y, g[3].z),
112+
y4m3(g[0].x, g[0].y, g[0].z), y4m3(g[1].x, g[1].y, g[1].z), y4m3(g[2].x, g[2].y, g[2].z), y4m3(g[3].x, g[3].y, g[3].z),
113+
y44(g[0].x, g[0].y, g[0].z), y44(g[1].x, g[1].y, g[1].z), y44(g[2].x, g[2].y, g[2].z), y44(g[3].x, g[3].y, g[3].z),
114+
y4m4(g[0].x, g[0].y, g[0].z), y4m4(g[1].x, g[1].y, g[1].z), y4m4(g[2].x, g[2].y, g[2].z), y4m4(g[3].x, g[3].y, g[3].z),
115+
0.000000000000000, 0.000000000000000, 0.935602579627389, 0.090028400200397,
116+
-0.452946651195697, -0.000000000000000, -0.000000000000000, -0.348678494661834,
117+
-0.000000000000000, -0.452946651195697, -0.000000000000000, -0.348678494661834,
118+
-0.000000000000000, 0.000000000000000, 0.000000000000000, -0.000000000000000,
119+
-0.000000000000000, -0.000000000000000, 0.000000000000000, -0.000000000000000,
120+
0.489238299435250, 0.000000000000000, -0.000000000000000, -0.376615818502422,
121+
0.000000000000000, -0.489238299435250, -0.000000000000000, 0.376615818502422,
122+
0.000000000000000, 0.000000000000000, 0.000000000000000, 0.532615198330370,
123+
0.000000000000000, 0.000000000000000, 0.000000000000000, -0.000000000000000,
124+
-0.656382056840170, -0.000000000000000, -0.000000000000000, -0.168427714314628,
125+
-0.000000000000000, -0.656382056840170, -0.000000000000000, -0.168427714314628,
126+
-0.317846011338142, -0.317846011338142, 1.017107236282055, 0.226023830284901,
127+
-0.000000000000000, -0.000000000000000, -0.000000000000000, 0.258942827786103,
128+
-0.000000000000000, -0.000000000000000, -0.000000000000000, 0.258942827786103,
129+
0.460602629757462, -0.460602629757462, 0.000000000000000, -0.000000000000000,
130+
0.000000000000000, 0.000000000000000, 0.000000000000000, -0.409424559784410,
131+
-0.000000000000000, -0.000000000000000, -0.000000000000000, 0.136474853261470,
132+
-0.000000000000000, 0.000000000000000, -0.000000000000000, -0.136474853261470,
133+
-0.504564900728724, -0.504564900728724, 0.000000000000000, -0.598002845308118,
134+
-0.000000000000000, -0.000000000000000, 0.000000000000000, 0.000000000000000,
135+
-0.000000000000000, -0.000000000000000, -0.000000000000000, 0.350610246256556,
136+
-0.000000000000000, -0.000000000000000, -0.000000000000000, 0.350610246256556,
137+
0.683184105191914, -0.683184105191914, 0.000000000000000, -0.000000000000000,
138+
0.000000000000000, 0.000000000000000, 0.000000000000000, -0.202424920056864,
139+
0.000000000000000, 0.000000000000000, 1.092548430592079, -0.350435072502801,
140+
0.451658037912587, 0.000000000000000, -0.000000000000000, 0.046358202625865,
141+
0.000000000000000, 0.451658037912587, -0.000000000000000, 0.046358202625865,
142+
0.000000000000000, -0.000000000000000, 0.000000000000000, 0.000000000000000,
143+
0.000000000000000, 0.000000000000000, 0.000000000000000, 0.492067081245654,
144+
-0.469376801586882, -0.000000000000000, -0.000000000000000, 0.187354445356332,
145+
-0.000000000000000, 0.469376801586882, -0.000000000000000, -0.187354445356332,
146+
0.000000000000000, 0.000000000000000, 0.000000000000000, 0.355076798886913,
147+
0.000000000000000, 0.000000000000000, 0.000000000000000, -0.000000000000000,
148+
0.518915578720260, 0.000000000000000, -0.000000000000000, -0.443845998608641,
149+
0.000000000000000, 0.518915578720260, -0.000000000000000, -0.443845998608641,
150+
0.000000000000000, -0.000000000000000, 0.000000000000000, 0.000000000000000,
151+
0.000000000000000, 0.000000000000000, 0.000000000000000, 0.452635881587108,
152+
-0.707162732524596, 0.000000000000000, -0.000000000000000, 0.120972027847095,
153+
-0.000000000000000, 0.707162732524596, -0.000000000000000, -0.120972027847095
154+
} ;
155+
156+
157+
}
158+
159+
void TearDown()
160+
{
161+
delete [] g;
162+
delete [] ref;
163+
delete [] rly;
164+
}
165+
};
166+
167+
TEST_F(YlmRealTest,YlmReal)
168+
{
169+
ModuleBase::YlmReal::Ylm_Real(nylm,ng,g,ylm);
170+
for(int i=0;i<nylm;++i)
171+
{
172+
for(int j=0;j<ng;++j)
173+
{
174+
EXPECT_NEAR(ylm(i,j),ref[i*ng+j],doublethreshold) << "i=" << i << " ,j=" << j;
175+
}
176+
}
177+
}
178+
179+
180+
TEST_F(YlmRealTest,YlmReal2)
181+
{
182+
ModuleBase::YlmReal::Ylm_Real2(nylm,ng,g,ylm);
183+
for(int i=0;i<nylm;++i)
184+
{
185+
for(int j=0;j<ng;++j)
186+
{
187+
EXPECT_NEAR(ylm(i,j),ref[i*ng+j],doublethreshold) << "i=" << i << " ,j=" << j;
188+
}
189+
}
190+
}
191+
192+
193+
TEST_F(YlmRealTest,rlylm)
194+
{
195+
for(int j=0;j<ng;++j)
196+
{
197+
ModuleBase::YlmReal::rlylm(lmax,g[j].x,g[j].y,g[j].z,rly);
198+
for(int i=0;i<nylm;++i)
199+
{
200+
EXPECT_NEAR(rly[i],ref[i*ng+j],doublethreshold) << "i=" << i << " ,j=" << j;
201+
}
202+
}
203+
}
204+

0 commit comments

Comments
 (0)