Skip to content

Commit 10c2079

Browse files
committed
Merge branch 'SDFT' of https://github.com/deepmodeling/abacus-develop into SDFT
2 parents b22d404 + aadec1b commit 10c2079

File tree

5 files changed

+257
-6
lines changed

5 files changed

+257
-6
lines changed

source/module_base/math_chebyshev.cpp

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,22 @@ Chebyshev<REAL>::~Chebyshev()
7676
delete [] coef_complex;
7777
}
7878

79+
template<typename REAL>
80+
void Chebyshev<REAL>::getpolyval(const REAL x, REAL* polyval, const int N)
81+
{
82+
polyval[0] = 1;
83+
polyval[1] = x;
84+
for(int i = 2; i < N; ++i)
85+
{
86+
polyval[i] = 2 * x * polyval[i-1] - polyval[i-2];
87+
}
88+
}
89+
template<typename REAL>
90+
REAL Chebyshev<REAL>::recurs(const REAL x, const REAL Tn, REAL const Tn_1)
91+
{
92+
return 2*x*Tn-Tn_1;
93+
}
94+
7995
template<typename REAL>
8096
REAL Chebyshev<REAL>:: ddot_real(
8197
const std::complex<REAL>* psi_L,

source/module_base/math_chebyshev.h

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ class FFTW;
2424
* T_{n+2}(x) = 2xT_{n+1}(x) - T_n(x)
2525
* II.
2626
* Any analytical function f(x) can be expanded by Chebyshev polynomial:
27-
* f(x) = \sum_{n=0}^{norder} C_n[f]*T_n(x) (|x| < 1),
27+
* f(x) = \sum_{n=0}^{norder-1} C_n[f]*T_n(x) (|x| < 1),
2828
* where C_n[f] = \frac{2-\delta_{n0}}{\pi} \int_0^\pi f(cos(\theta))cos(n\theta) d\theta
2929
* Here C_n can be calculate with FFT.
3030
* III.
@@ -40,13 +40,13 @@ class FFTW;
4040
* USAGE:
4141
* Chebyshev che(10); // constructe a chebyshev expansion of 10 orders (n=0,1,...,9)
4242
* 1. che.calcoef_real(&a, &A::cos) // calculate C_n[f], where f is a.cos
43-
* for(inti=0;i<10;++i) cout<<che.coef_real[i]<<endl; //Then we print C_n[f]
43+
* for(int i=0;i<10;++i) cout<<che.coef_real[i]<<endl; //Then we print C_n[f]
4444
*
4545
* che.calcoef_complex(&b, &B::expi) // calculate C_n[g], where g is b.expi
46-
* for(inti=0;i<10;++i) cout<<che.coef_complex[i]<<endl; //Then we print C_n[g]
46+
* for(int i=0;i<10;++i) cout<<che.coef_complex[i]<<endl; //Then we print C_n[g]
4747
*
4848
* che.calcoef_pair(&c, &C::cos, &C::sin) // calculate C_n[g], where g is (c.cos, c.sin)
49-
* for(inti=0;i<10;++i) cout<<che.coef_complex[i]<<endl; //Then we print C_n[g]
49+
* for(int i=0;i<10;++i) cout<<che.coef_complex[i]<<endl; //Then we print C_n[g]
5050
*
5151
* 2. che.calcoef_real(&occ, &Occupy::fd)
5252
* che.calfinalvec_real(&hamilt, &Hamilt::hpsi, psi_in, psi_out, npw);
@@ -59,7 +59,15 @@ class FFTW;
5959
* 3. che.tracepolyA(&hamilt, &Hamilt::hpsi, psi_in, npw, npwx, nbands)
6060
* //calculate \sum_i^{nbands} <psi_i|T_n(H)|psi_i>
6161
*
62-
* 4. che.recurs_complex(&hamilt, &Hamilt::hpsi, vp1, v, vm1, npw)
62+
* 4. che.calcoef_complex(&fun, &toolfunc::expi); //calculate C_n[exp(ix)]
63+
* che.getpolyval(PI/4, T, norder); //get T_n(pi/4)
64+
* std::complex<double> sum(0,0);
65+
* for(int i = 0; i < norder ; ++i)
66+
* {
67+
* sum += che.coef_complex[i]*T[i]; //sum = exp(i*pi/4) = \sum_n C_n[exp(ix)]*T_n(pi/4)
68+
* }
69+
*
70+
* 5. che.recurs_complex(&hamilt, &Hamilt::hpsi, vp1, v, vm1, npw)
6371
* //calculate vp1: |vp1> = 2 H|v> - |vm1>;
6472
*
6573
*/
@@ -121,6 +129,8 @@ class Chebyshev
121129
T *ptr, void (T::*funA)(std::complex<REAL> *in, std::complex<REAL> *out, const int),
122130
std::complex<REAL> *wavein,
123131
const int N, const int LDA = 1, const int m = 1);
132+
//get T_n(x)
133+
void getpolyval(REAL x, REAL* polyval, const int N);
124134

125135
// IV.
126136
// recurs fomula: v_{n+1} = 2Av_n - v_{n-1}
@@ -141,6 +151,8 @@ class Chebyshev
141151
std::complex<REAL>* arrayn, //v_n
142152
std::complex<REAL>* arrayn_1, //v_{n-1}
143153
const int N, const int LDA = 1, const int m = 1);
154+
//return 2xTn-Tn_1
155+
inline REAL recurs(const REAL x, const REAL Tn, const REAL Tn_1);
144156

145157
// V.
146158
// auxiliary function

source/module_base/math_chebyshev_def.h

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -429,9 +429,13 @@ bool Chebyshev<REAL>::checkconverge(
429429
(ptr->*funA)(arrayn_1, arrayn,1);
430430
REAL sum1,sum2;
431431
REAL t;
432-
432+
#ifdef __MPI
433433
sum1=ModuleBase::GlobalFunc::ddot_real(N,arrayn_1,arrayn_1);
434434
sum2=ModuleBase::GlobalFunc::ddot_real(N,arrayn_1,arrayn);
435+
#else
436+
sum1=this->ddot_real(arrayn_1,arrayn_1,N);
437+
sum2=this->ddot_real(arrayn_1,arrayn,N);
438+
#endif
435439
t = sum2 / sum1 * (tmax - tmin) / 2 + (tmax + tmin) / 2;
436440
if(t < tmin || tmin == 0)
437441
{
@@ -447,8 +451,13 @@ bool Chebyshev<REAL>::checkconverge(
447451
for(int ior = 2; ior < norder; ++ior)
448452
{
449453
(ptr->*funA)(arrayn,arraynp1,1);
454+
#ifdef __MPI
450455
sum1=ModuleBase::GlobalFunc::ddot_real(N,arrayn,arrayn);
451456
sum2=ModuleBase::GlobalFunc::ddot_real(N,arrayn,arraynp1);
457+
#else
458+
sum1=this->ddot_real(arrayn,arrayn,N);
459+
sum2=this->ddot_real(arrayn,arraynp1,N);
460+
#endif
452461
t = sum2/sum1 * (tmax - tmin) / 2 + (tmax + tmin) / 2;
453462
if(t < tmin)
454463
{

source/module_base/test/CMakeLists.txt

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -108,3 +108,9 @@ AddTest(
108108
TARGET base_container
109109
SOURCES container_operator_test.cpp ../container_operator.h
110110
)
111+
112+
AddTest(
113+
TARGET base_math_chebyshev
114+
LIBS ${math_libs}
115+
SOURCES math_chebyshev_test.cpp ../math_chebyshev.cpp ../global_function.cpp ../tool_quit.cpp ../global_variable.cpp ../timer.cpp ../global_file.cpp ../memory.cpp
116+
)
Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
#include "../math_chebyshev.h"
2+
#include"gtest/gtest.h"
3+
/************************************************
4+
* unit test of class Chebyshev
5+
***********************************************/
6+
7+
/**
8+
* - Tested Functions:
9+
* - calcoef_real
10+
* - calcoef_complex
11+
* - calcoef_pair
12+
* - calfinalvec_real
13+
* - calfinalvec_complex
14+
* - tracepolyA
15+
* - checkconverge
16+
*
17+
*
18+
*/
19+
class toolfunc{
20+
public:
21+
double x7(double x){
22+
return pow(x,7);
23+
}
24+
double x6(double x){
25+
return pow(x,6);
26+
}
27+
double expr(double x){
28+
return exp(x);
29+
}
30+
std::complex<double> expi(std::complex<double> x){
31+
const std::complex<double> j(0.0,1.0);
32+
return exp(j*x);
33+
}
34+
std::complex<double> expi2(std::complex<double> x){
35+
const std::complex<double> j(0.0,1.0);
36+
const double PI = 3.14159265358979323846;
37+
return exp(j*PI/2.0*x);
38+
}
39+
//Pauli matrix: [0,-i;i,0]
40+
void sigma_y(std::complex<double>* spin_in, std::complex<double>* spin_out, const int m = 1){
41+
const std::complex<double> j(0.0,1.0);
42+
for(int i = 0 ; i < m ; ++i)
43+
{
44+
spin_out[2*i] = -j*spin_in[2*i+1];
45+
spin_out[2*i+1] = j*spin_in[2*i];
46+
}
47+
}
48+
};
49+
class MathChebyshevTest : public testing::Test
50+
{
51+
protected:
52+
ModuleBase::Chebyshev<double> *p_chetest;
53+
toolfunc fun;
54+
};
55+
56+
57+
TEST_F(MathChebyshevTest,calcoef_real)
58+
{
59+
p_chetest = new ModuleBase::Chebyshev<double>(10);
60+
//x^6 = 1/32*( 10T_0 + 15T_2 + 6T_4 + T_6 )
61+
//x^7 = 1/64*( 35T_1 + 21T_3 + 7T_5 + T_7 )
62+
const double x6ref[10] = {10,0,15,0,6,0,1,0,0,0};
63+
const double x7ref[10] = {0,35,0,21,0,7,0,1,0,0};
64+
p_chetest->calcoef_real(&fun,&toolfunc::x6);
65+
for(int i = 0;i < 10 ; ++i)
66+
{
67+
EXPECT_NEAR(p_chetest->coef_real[i]*32.0, x6ref[i] ,1.e-8);
68+
}
69+
p_chetest->calcoef_real(&fun,&toolfunc::x7);
70+
for(int i = 0;i < 10 ; ++i)
71+
{
72+
EXPECT_NEAR(p_chetest->coef_real[i]*64.0, x7ref[i] ,1.e-8);
73+
}
74+
delete p_chetest;
75+
}
76+
77+
TEST_F(MathChebyshevTest,calcoef_pair)
78+
{
79+
p_chetest = new ModuleBase::Chebyshev<double>(10);
80+
//x^6 = 1/32*( 10T_0 + 15T_2 + 6T_4 + T_6 )
81+
//x^7 = 1/64*( 35T_1 + 21T_3 + 7T_5 + T_7 )
82+
const double x6ref[10] = {10,0,15,0,6,0,1,0,0,0};
83+
const double x7ref[10] = {0,35,0,21,0,7,0,1,0,0};
84+
p_chetest->calcoef_pair(&fun, &toolfunc::x6, &toolfunc::x7);
85+
for(int i = 0;i < 10 ; ++i)
86+
{
87+
EXPECT_NEAR(p_chetest->coef_complex[i].real()*32.0, x6ref[i] ,1.e-8);
88+
EXPECT_NEAR(p_chetest->coef_complex[i].imag()*64.0, x7ref[i] ,1.e-8);
89+
}
90+
delete p_chetest;
91+
}
92+
93+
TEST_F(MathChebyshevTest,calcoef_complex)
94+
{
95+
const int norder = 100;
96+
const double PI = 3.14159265358979323846;
97+
p_chetest = new ModuleBase::Chebyshev<double>(norder);
98+
double *T = new double [norder];
99+
//check exp(i\pi/4) = \sum_n C_n[exp(ix)]T_n(\pi/4) = sqrt(2)/2*(1, i)
100+
p_chetest->calcoef_complex(&fun, &toolfunc::expi);
101+
p_chetest->getpolyval(PI/4, T, norder);
102+
std::complex<double> sum(0,0);
103+
for(int i = 0; i < norder ; ++i)
104+
{
105+
sum += p_chetest->coef_complex[i]*T[i];
106+
}
107+
EXPECT_NEAR(sum.real(), sqrt(2)/2 ,1.e-8);
108+
EXPECT_NEAR(sum.imag(), sqrt(2)/2 ,1.e-8);
109+
delete []T;
110+
delete p_chetest;
111+
}
112+
113+
TEST_F(MathChebyshevTest,calfinalvec_real)
114+
{
115+
const int norder = 100;
116+
const double E = 2.718281828459046;
117+
p_chetest = new ModuleBase::Chebyshev<double>(norder);
118+
// 1 [ 1/e+e -i(e-1/e) ]
119+
// exp(\sigma_y)= - [ ], where \sigma_y = [0, -i; i, 0]
120+
// 2 [ i(e-1/e) 1/e+e ]
121+
std::complex<double> *v = new std::complex<double> [4];
122+
std::complex<double> *vout = new std::complex<double> [4];
123+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
124+
125+
p_chetest->calcoef_real(&fun,&toolfunc::expr);
126+
p_chetest->calfinalvec_real(&fun, &toolfunc::sigma_y, v, vout, 2,2,2);
127+
EXPECT_NEAR(vout[0].real(), 0.5*(E+1/E) ,1.e-8);
128+
EXPECT_NEAR(vout[0].imag(), 0 ,1.e-8);
129+
EXPECT_NEAR(vout[1].real(), 0 ,1.e-8);
130+
EXPECT_NEAR(vout[1].imag(), 0.5*(E-1/E) ,1.e-8);
131+
EXPECT_NEAR(vout[2].real(), 0 ,1.e-8);
132+
EXPECT_NEAR(vout[2].imag(),-0.5*(E-1/E) ,1.e-8);
133+
EXPECT_NEAR(vout[3].real(), 0.5*(E+1/E) ,1.e-8);
134+
EXPECT_NEAR(vout[3].imag(), 0 ,1.e-8);
135+
136+
delete[] v;
137+
delete[] vout;
138+
delete p_chetest;
139+
}
140+
141+
TEST_F(MathChebyshevTest,calfinalvec_complex)
142+
{
143+
const int norder = 100;
144+
const double E = 2.718281828459046;
145+
p_chetest = new ModuleBase::Chebyshev<double>(norder);
146+
// [ 0 1 ]
147+
// exp(i pi/2*\sigma_y)= [ ], where \sigma_y = [0, -i; i, 0]
148+
// [ -1 0 ]
149+
std::complex<double> *v = new std::complex<double> [4];
150+
std::complex<double> *vout = new std::complex<double> [4];
151+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
152+
153+
p_chetest->calcoef_complex(&fun,&toolfunc::expi2);
154+
p_chetest->calfinalvec_complex(&fun, &toolfunc::sigma_y, v, vout, 2,2,2);
155+
EXPECT_NEAR(vout[0].real(), 0 ,1.e-8);
156+
EXPECT_NEAR(vout[0].imag(), 0 ,1.e-8);
157+
EXPECT_NEAR(vout[1].real(),-1 ,1.e-8);
158+
EXPECT_NEAR(vout[1].imag(), 0 ,1.e-8);
159+
EXPECT_NEAR(vout[2].real(), 1 ,1.e-8);
160+
EXPECT_NEAR(vout[2].imag(), 0 ,1.e-8);
161+
EXPECT_NEAR(vout[3].real(), 0 ,1.e-8);
162+
EXPECT_NEAR(vout[3].imag(), 0 ,1.e-8);
163+
164+
delete[] v;
165+
delete[] vout;
166+
delete p_chetest;
167+
}
168+
169+
TEST_F(MathChebyshevTest,tracepolyA)
170+
{
171+
const int norder = 100;
172+
p_chetest = new ModuleBase::Chebyshev<double>(norder);
173+
174+
std::complex<double> *v = new std::complex<double> [4];
175+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
176+
177+
p_chetest->tracepolyA(&fun, &toolfunc::sigma_y, v, 2,2,2);
178+
//Trace: even function: 2 ; odd function 0.
179+
for(int i = 0 ; i < norder ; ++i)
180+
{
181+
if(i%2==0) EXPECT_NEAR(p_chetest->polytrace[i], 2 ,1.e-8);
182+
else EXPECT_NEAR(p_chetest->polytrace[i], 0 ,1.e-8);
183+
}
184+
185+
delete[] v;
186+
delete p_chetest;
187+
}
188+
189+
TEST_F(MathChebyshevTest,checkconverge)
190+
{
191+
const int norder = 100;
192+
p_chetest = new ModuleBase::Chebyshev<double>(norder);
193+
194+
std::complex<double> *v = new std::complex<double> [4];
195+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
196+
double tmin = -1.1;
197+
double tmax = 1.1;
198+
bool converge;
199+
converge = p_chetest->checkconverge(&fun, &toolfunc::sigma_y, v, 2, tmax, tmin, 0.2);
200+
EXPECT_TRUE(converge);
201+
converge = p_chetest->checkconverge(&fun, &toolfunc::sigma_y, v+2, 2, tmax, tmin, 0.2);
202+
EXPECT_TRUE(converge);
203+
EXPECT_NEAR(tmin, -1.1, 1e-8);
204+
EXPECT_NEAR(tmax, 1.1, 1e-8);
205+
206+
delete[] v;
207+
delete p_chetest;
208+
}

0 commit comments

Comments
 (0)