Skip to content

Commit bbf41aa

Browse files
YuLiu98YuLiu98
authored
Refactor: replace 3.14 using PI (#6085)
* Refactor: replace 3.14 using PI * Test: update result.ref * Test: update result.ref
1 parent 60a408a commit bbf41aa

File tree

14 files changed

+83
-166
lines changed

14 files changed

+83
-166
lines changed

‎source/module_basis/module_ao/ORB_nonlocal_lm.cpp‎

Lines changed: 13 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,18 @@
1-
#include <cassert>
2-
#include <iostream>
3-
#include "module_parameter/parameter.h"
4-
#include <sstream>
5-
#include <cmath>
61
#include "ORB_nonlocal_lm.h"
7-
#include "module_base/math_integral.h"
2+
3+
#include "module_base/constants.h"
84
#include "module_base/global_function.h"
9-
#include "module_base/mathzone.h" /// use Polynomial_Interpolation_xy, Spherical_Bessel
5+
#include "module_base/math_integral.h"
6+
#include "module_base/math_polyint.h"
7+
#include "module_base/math_sphbes.h"
8+
#include "module_base/mathzone.h" /// use Polynomial_Interpolation_xy, Spherical_Bessel
109
#include "module_base/mathzone_add1.h" /// use SplineD2
11-
#include "module_base/math_sphbes.h" // mohan add 2021-05-06
12-
#include "module_base/math_polyint.h" // mohan add 2021-05-06
10+
#include "module_parameter/parameter.h"
11+
12+
#include <cassert>
13+
#include <cmath>
14+
#include <iostream>
15+
#include <sstream>
1316

1417
Numerical_Nonlocal_Lm::Numerical_Nonlocal_Lm()
1518
{
@@ -178,84 +181,13 @@ void Numerical_Nonlocal_Lm::set_NL_proj(
178181
return;
179182
}
180183

181-
/*
182-
// extra_uniform currently does not work, because
183-
// beta[ir] = this->beta_r[ir]/r_radial[ir];
184-
// does not work for ir==0, and this formula is not stable for small r.
185-
//
186-
// If one really needs beta (in what circumstance?), one should do
187-
// extrapolation to reach r=0.
188-
189-
void Numerical_Nonlocal_Lm::extra_uniform(const double &dr_uniform_in)
190-
{
191-
assert(dr_uniform_in>0.0);
192-
this->dr_uniform = dr_uniform_in;
193-
this->nr_uniform = static_cast<int>(rcut/dr_uniform) + 10;
194-
195-
// std::cout << " nr_uniform = " << nr_uniform << std::endl;
196-
197-
delete[] this->beta_uniform;
198-
this->beta_uniform = new double[nr_uniform];
199-
ModuleBase::GlobalFunc::ZEROS (this->beta_uniform, nr_uniform);
200-
201-
// mohan fix bug 2011-04-14
202-
// the beta_r is beta*r.
203-
// and beta is beta!!!
204-
double* beta = new double[nr];
205-
ModuleBase::GlobalFunc::ZEROS(beta,nr);
206-
for(int ir=0; ir<nr; ir++)
207-
{
208-
assert(r_radial[ir]>0.0);
209-
beta[ir] = this->beta_r[ir]/r_radial[ir];
210-
}
211-
212-
for (int ir = 0; ir < this->nr_uniform; ir++)
213-
{
214-
double rnew = ir * dr_uniform;
215-
this->beta_uniform[ir] = ModuleBase::PolyInt::Polynomial_Interpolation_xy(this->r_radial, beta, this->nr, rnew);
216-
}
217-
delete[] beta;
218-
219-
delete[] this->dbeta_uniform;
220-
this->dbeta_uniform = new double[nr_uniform];
221-
ModuleBase::GlobalFunc::ZEROS(this->dbeta_uniform, nr_uniform);
222-
223-
double* y2 = new double[nr];
224-
ModuleBase::Mathzone_Add1::SplineD2 (r_radial, beta_r, nr, 0.0, 0.0, y2);
225-
226-
double* rad = new double[nr_uniform];
227-
for (int ir = 0; ir < nr_uniform; ir++)
228-
{
229-
rad[ir] = ir*dr_uniform;
230-
}
231-
232-
double* tmp = new double[nr_uniform];
233-
double* tmp1 = new double[nr_uniform];
234-
ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(r_radial, beta_r, y2,
235-
nr, rad, nr_uniform, tmp, dbeta_uniform);
236-
237-
for(int ir= 0 ; ir<nr_uniform; ir++)
238-
{
239-
//assert(dbeta_uniform[ir]==beta_uniform[ir]);
240-
}
241-
242-
delete [] y2;
243-
delete [] rad;
244-
delete [] tmp;
245-
delete [] tmp1;
246-
return;
247-
}
248-
*/
249-
250184
void Numerical_Nonlocal_Lm::get_kradial()
251185
{
252186
//ModuleBase::TITLE("Numerical_Nonlocal_Lm","get_kradial");
253187
double *jl = new double[nr];
254188
double *integrated_func = new double[nr];
255189

256-
// mohan add constant of pi, 2021-04-26
257-
const double pi = 3.14159265358979323846;
258-
const double pref = sqrt( 2.0 / pi );
190+
const double pref = sqrt(2.0 / ModuleBase::PI);
259191

260192
for (int ik = 0; ik < nk; ik++)
261193
{

‎source/module_cell/module_neighbor/sltk_util.h‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,6 @@
77

88
static struct { template<typename T> operator T*() const { return static_cast<T*>(0); } } NullPtr;
99

10-
static const double Pi = 3.1415926535897932384626433832795;
11-
12-
1310
/*** Function ***/
1411
template<typename exception>
1512
static inline void affirm(const bool b)

‎source/module_cell/module_paw/paw_cell.cpp‎

Lines changed: 7 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,9 @@
11
#include "paw_cell.h"
2+
3+
#include "module_base/constants.h"
4+
#include "module_base/tool_quit.h"
25
#include "module_base/tool_title.h"
36
#include "module_parameter/parameter.h"
4-
#include "module_base/tool_quit.h"
57
#ifdef __MPI
68
#include "module_base/parallel_reduce.h"
79
#endif
@@ -176,9 +178,6 @@ void Paw_Cell::set_paw_k(
176178
{
177179
ModuleBase::TITLE("Paw_Element","set_paw_k");
178180

179-
const double pi = 3.141592653589793238462643383279502884197;
180-
const double twopi = 2.0 * pi;
181-
182181
this -> npw = npw_in;
183182
this -> npwx = npwx_in;
184183

@@ -190,7 +189,7 @@ void Paw_Cell::set_paw_k(
190189
{
191190
arg += atom_coord[iat][i] * kpt[i];
192191
}
193-
arg *= twopi;
192+
arg *= ModuleBase::TWO_PI;
194193
const std::complex<double> kphase = std::complex<double>(cos(arg), -sin(arg));
195194

196195
struc_fact[iat].resize(npw);
@@ -326,17 +325,14 @@ std::vector<double> Paw_Cell::calc_ylm(const int lmax, const double * r)
326325
std::vector<std::complex<double>> phase;
327326
phase.resize(lmax+1);
328327

329-
const double pi = 3.141592653589793238462643383279502884197;
330-
const double fourpi = 4.0 * pi;
331-
332328
// set zero
333329
for(int i = 0; i < size_ylm; i++)
334330
{
335331
ylm[i] = 0;
336332
}
337333

338334
// l = 0
339-
ylm[0] = 1.0/std::sqrt(fourpi);
335+
ylm[0] = 1.0 / std::sqrt(ModuleBase::FOUR_PI);
340336

341337
if(rr>tol)
342338
{
@@ -361,8 +357,8 @@ std::vector<double> Paw_Cell::calc_ylm(const int lmax, const double * r)
361357
{
362358
const int l0 = l*l + l;
363359
double fact = 1.0/double(l*(l+1));
364-
const double ylmcst = std::sqrt(double(2*l+1)/fourpi);
365-
360+
const double ylmcst = std::sqrt(double(2 * l + 1) / ModuleBase::FOUR_PI);
361+
366362
// m = 0
367363
ylm[l0] = ylmcst*Paw_Cell::ass_leg_pol(l,0,ctheta);
368364

‎source/module_cell/module_paw/paw_sphbes.cpp‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,12 @@
1-
#include "paw_element.h"
2-
#include "module_base/tool_title.h"
1+
#include "module_base/constants.h"
32
#include "module_base/tool_quit.h"
3+
#include "module_base/tool_title.h"
4+
#include "paw_element.h"
45

56
double Paw_Element::get_ptilde(const int istate_in, const double q_in, const double omega)
67
{
78
// multiply by a factor 4pi / sqrt(omega)
8-
const double factor = 4.0 * 3.141592653589793238462643383279502884197 / std::sqrt(omega);
9+
const double factor = ModuleBase::FOUR_PI / std::sqrt(omega);
910

1011
return this->splint(qgrid, ptilde_q[istate_in], d2ptilde_q[istate_in], q_in) * factor;
1112
}
@@ -18,9 +19,6 @@ void Paw_Element::transform_ptilde()
1819
double dq = 0.01;
1920
int nq = int( (std::sqrt(ecutwfc) / dq + 4) * cell_factor );
2021

21-
const double pi = 3.141592653589793238462643383279502884197;
22-
const double twopi = 2.0 * pi;
23-
2422
std::vector<double> integrand;
2523
integrand.resize(nr);
2624

@@ -49,7 +47,7 @@ void Paw_Element::transform_ptilde()
4947
{
5048
for(int ir = 0; ir < nr; ir ++)
5149
{
52-
integrand[ir] = twopi * ptilde_r[istate][ir] * std::pow(rr[ir],3);
50+
integrand[ir] = ModuleBase::TWO_PI * ptilde_r[istate][ir] * std::pow(rr[ir], 3);
5351
}
5452
yp1 = this -> simpson_integration(integrand) / 3.0;
5553
}
@@ -60,7 +58,7 @@ void Paw_Element::transform_ptilde()
6058
double x = rr[ir] * double(nq - 1) * dq;
6159

6260
this -> spherical_bessel_function(l,x,bes,besp,1);
63-
integrand[ir] = twopi * besp * ptilde_r[istate][ir] * std::pow(rr[ir],3);
61+
integrand[ir] = ModuleBase::TWO_PI * besp * ptilde_r[istate][ir] * std::pow(rr[ir], 3);
6462
}
6563

6664
ypn = this -> simpson_integration(integrand);

‎source/module_elecstate/module_charge/charge_init.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -196,8 +196,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
196196
{
197197
std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
198198
}
199-
const double pi = 3.141592653589790;
200-
const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
199+
const double fact = (3.0 / 5.0) * pow(3.0 * ModuleBase::PI * ModuleBase::PI, 2.0 / 3.0);
201200
for (int is = 0; is < PARAM.inp.nspin; ++is)
202201
{
203202
for (int ir = 0; ir < this->rhopw->nrxx; ++ir)

‎source/module_hamilt_general/module_xc/xc_funct_exch_gga.cpp‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -96,8 +96,7 @@ double &sx, double &v1x, double &v2x)
9696

9797
// numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
9898
const double third = 1.0 / 3.0;
99-
const double pi = 3.14159265358979323846;
100-
const double c1 = 0.750 / pi;
99+
const double c1 = 0.750 / ModuleBase::PI;
101100
const double c2 = 3.0936677262801360;
102101
const double c5 = 4.0 * third;
103102
// parameters of the functional
@@ -196,11 +195,10 @@ void XC_Functional::wcx(const double &rho,const double &grho, double &sx, double
196195
// exchange energy gradient part
197196
double dxunif, dfx, f1, f2, f3, dfx1, x1, x2, x3, dxds1, dxds2, dxds3;
198197
// numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
199-
double third, c1, c2, c5, teneightyone; // c6
200-
const double pi = 3.14159265358979323846;
198+
double third, c1, c2, c5, teneightyone; // c6
201199

202200
third = 1.0/3.0;
203-
c1 = 0.75/pi;
201+
c1 = 0.75 / ModuleBase::PI;
204202
c2 = 3.093667726280136;
205203
c5 = 4.0 * third;
206204
teneightyone = 0.123456790123;

‎source/module_hamilt_general/module_xc/xc_funct_exch_lda.cpp‎

Lines changed: 8 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -35,15 +35,14 @@ void XC_Functional::slater1(const double &rs, double &ex, double &vx)
3535
// Slater exchange with alpha=2/3 and Relativistic exchange
3636
void XC_Functional::slater_rxc(const double &rs, double &ex, double &vx)
3737
{
38-
static const double pi = 3.14159265358979;
3938
const double trd = 1.0 / 3.0;
4039
//const double ftrd = 4.0 / 3.0;
4140
//const double tftm = pow(2.0, ftrd) - 2.0;
42-
const double a0 = pow((4.0 / (9.0 * pi)), trd);
41+
const double a0 = pow((4.0 / (9.0 * ModuleBase::PI)), trd);
4342
// X-alpha parameter:
4443
const double alp = 2 * trd;
4544

46-
double vxp = -3 * alp / (2 * pi * a0 * rs);
45+
double vxp = -3 * alp / (2 * ModuleBase::PI * a0 * rs);
4746
double exp = 3 * vxp / 4;
4847
const double beta = 0.014 / rs;
4948
const double sb = sqrt(1 + beta * beta);
@@ -53,7 +52,7 @@ void XC_Functional::slater_rxc(const double &rs, double &ex, double &vx)
5352
exp = exp * (1.0 - 1.5 * x * x);
5453
vx = vxp;
5554
ex = exp;
56-
return;
55+
return;
5756
}
5857

5958
// Slater exchange with alpha=2/3, spin-polarized case
@@ -107,20 +106,18 @@ void XC_Functional::slater_rxc_spin( const double &rho, const double &z,
107106
return;
108107
}
109108

110-
const double pi = 3.141592653589790;
111-
112-
const double trd = 1.0 / 3.0;
109+
const double trd = 1.0 / 3.0;
113110
const double ftrd = 4.0 / 3.0;
114111
double tftm = pow(2.0, ftrd) - 2;
115-
double a0 = pow((4 / (9 * pi)), trd);
112+
double a0 = pow((4 / (9 * ModuleBase::PI)), trd);
116113

117-
double alp = 2 * trd;
114+
double alp = 2 * trd;
118115

119116
double fz = (pow((1 + z), ftrd) + pow((1 - z), ftrd) - 2) / tftm;
120117
double fzp = ftrd * (pow((1 + z), trd) - pow((1 - z), trd)) / tftm;
121118

122-
double rs = pow((3 / (4 * pi * rho)), trd);
123-
double vxp = -3 * alp / (2 * pi * a0 * rs);
119+
double rs = pow((3 / (4 * ModuleBase::PI * rho)), trd);
120+
double vxp = -3 * alp / (2 * ModuleBase::PI * a0 * rs);
124121
double exp = 3 * vxp / 4;
125122

126123
double beta = 0.014 / rs;

‎source/module_hamilt_general/module_xc/xc_funct_hcth.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,6 @@ void XC_Functional::hcth(const double rho, const double grho, double &sx, double
2020
// real(kind=DP) :: rho, grho, sx, v1x, v2x
2121

2222
// parameter :
23-
const double pi = 3.14159265358979323846;
2423
double o3 = 1.00 / 3.00;
2524
double o34 = 4.00 / 3.00;
2625
double fr83 = 8.0 / 3.0;
@@ -33,7 +32,7 @@ void XC_Functional::hcth(const double rho, const double grho, double &sx, double
3332
dgaa_drho, dgab_drho, dgx_drho, dgaa_dgr, dgab_dgr, dgx_dgr;
3433

3534
r3q2 = std::pow(2.0, (-o3));
36-
r3pi = std::pow((3.0 / pi), o3);
35+
r3pi = std::pow((3.0 / ModuleBase::PI), o3);
3736
//.....coefficients for pwf correlation......................................
3837
cg0[1] = 0.0310910;
3938
cg0[2] = 0.2137000;

‎source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp‎

Lines changed: 1 addition & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -176,19 +176,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
176176

177177
double *aux1 = new double[msh];
178178

179-
// for tests
180-
/*
181-
for(ir=0; ir<msh; ir++)
182-
{
183-
aux[ir] = r[ir] * zp_in * e2 / ucell.omega;
184-
}
185-
ModuleBase::Integral::Simpson_Integral(msh, aux, rab, vloc_1d[0] );
186-
vloc_1d[0] *= 4*3.1415926;
187-
std::cout << " vloc_1d[0]=" << vloc_1d[0]/rho_basis->npw << std::endl;
188-
std::cout << " vloc_1d[0]=" << vloc_1d[0]/rho_basis->nxyz << std::endl;
189-
*/
190-
191-
// (1)
179+
// (1)
192180
if(rho_basis->gg_uniq[0] < 1.0e-8)
193181
{
194182
double *aux = new double[msh];

0 commit comments

Comments
 (0)