Refactor: replace 3.14 using PI

source/module_basis/module_ao/ORB_nonlocal_lm.cpp

@@ -1,15 +1,18 @@
-#include <cassert>
-#include <iostream>
-#include "module_parameter/parameter.h"
-#include <sstream>
-#include <cmath>
 #include "ORB_nonlocal_lm.h"
-#include "module_base/math_integral.h"
+
+#include "module_base/constants.h"
 #include "module_base/global_function.h"
-#include "module_base/mathzone.h" /// use Polynomial_Interpolation_xy, Spherical_Bessel
+#include "module_base/math_integral.h"
+#include "module_base/math_polyint.h"
+#include "module_base/math_sphbes.h"
+#include "module_base/mathzone.h"      /// use Polynomial_Interpolation_xy, Spherical_Bessel
 #include "module_base/mathzone_add1.h" /// use SplineD2
-#include "module_base/math_sphbes.h" // mohan add 2021-05-06
-#include "module_base/math_polyint.h" // mohan add 2021-05-06
+#include "module_parameter/parameter.h"
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <sstream>
 
 Numerical_Nonlocal_Lm::Numerical_Nonlocal_Lm()
 {
@@ -178,84 +181,13 @@ void Numerical_Nonlocal_Lm::set_NL_proj(
 	return;
 }
 
-/*
-// extra_uniform currently does not work, because
-// beta[ir] = this->beta_r[ir]/r_radial[ir];
-// does not work for ir==0, and this formula is not stable for small r.
-//
-// If one really needs beta (in what circumstance?), one should do
-// extrapolation to reach r=0.
-
-void Numerical_Nonlocal_Lm::extra_uniform(const double &dr_uniform_in)
-{
-	assert(dr_uniform_in>0.0);
-	this->dr_uniform = dr_uniform_in;
-	this->nr_uniform = static_cast<int>(rcut/dr_uniform) + 10;
-
-//	std::cout << " nr_uniform = " << nr_uniform << std::endl;
-
-	delete[] this->beta_uniform;
-	this->beta_uniform = new double[nr_uniform];
-	ModuleBase::GlobalFunc::ZEROS (this->beta_uniform, nr_uniform);
-
-	// mohan fix bug 2011-04-14
-	// the beta_r is beta*r.
-	// and beta is beta!!!
-	double* beta = new double[nr];
-	ModuleBase::GlobalFunc::ZEROS(beta,nr);
-	for(int ir=0; ir<nr; ir++)
-	{
-		assert(r_radial[ir]>0.0);
-		beta[ir] = this->beta_r[ir]/r_radial[ir];
-	}
-
-	for (int ir = 0; ir < this->nr_uniform; ir++)
-	{
-		double rnew = ir * dr_uniform;
-		this->beta_uniform[ir] = ModuleBase::PolyInt::Polynomial_Interpolation_xy(this->r_radial, beta, this->nr, rnew); 
-    }
-	delete[] beta;
-
-	delete[] this->dbeta_uniform;
-	this->dbeta_uniform = new double[nr_uniform];
-	ModuleBase::GlobalFunc::ZEROS(this->dbeta_uniform, nr_uniform);
-
-	double* y2 = new double[nr];
-	ModuleBase::Mathzone_Add1::SplineD2 (r_radial, beta_r, nr, 0.0, 0.0, y2);
-
-	double* rad = new double[nr_uniform];
-	for (int ir = 0; ir < nr_uniform; ir++)
-	{
-		rad[ir] = ir*dr_uniform;
-	}
-
-	double* tmp = new double[nr_uniform];
-	double* tmp1 = new double[nr_uniform];
-	ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(r_radial, beta_r, y2, 
-					nr, rad, nr_uniform, tmp, dbeta_uniform);
-
-	for(int ir= 0 ; ir<nr_uniform; ir++)
-	{
-		//assert(dbeta_uniform[ir]==beta_uniform[ir]);
-	}
-
-	delete [] y2;
-	delete [] rad;	
-	delete [] tmp;
-	delete [] tmp1;
-	return;
-}
-*/
-
 void Numerical_Nonlocal_Lm::get_kradial()
 {
     //ModuleBase::TITLE("Numerical_Nonlocal_Lm","get_kradial");
     double *jl = new double[nr];
     double *integrated_func = new double[nr];
 
-	// mohan add constant of pi, 2021-04-26
-	const double pi = 3.14159265358979323846;
-    const double pref = sqrt( 2.0 / pi );
+    const double pref = sqrt(2.0 / ModuleBase::PI);
 
     for (int ik = 0; ik < nk; ik++)
     {

source/module_cell/module_neighbor/sltk_util.h

@@ -7,9 +7,6 @@
 
 static struct { template<typename T> operator T*() const { return static_cast<T*>(0); } } NullPtr;
 
-static const double Pi = 3.1415926535897932384626433832795;
-
-
 /*** Function ***/
 template<typename exception>
 static inline void affirm(const bool b)

source/module_cell/module_paw/paw_cell.cpp

@@ -1,7 +1,9 @@
 #include "paw_cell.h"
+
+#include "module_base/constants.h"
+#include "module_base/tool_quit.h"
 #include "module_base/tool_title.h"
 #include "module_parameter/parameter.h"
-#include "module_base/tool_quit.h"
 #ifdef __MPI
 #include "module_base/parallel_reduce.h"
 #endif
@@ -176,9 +178,6 @@ void Paw_Cell::set_paw_k(
 {
     ModuleBase::TITLE("Paw_Element","set_paw_k");
 
-    const double pi = 3.141592653589793238462643383279502884197;
-    const double twopi = 2.0 * pi;
-
     this -> npw = npw_in;
     this -> npwx = npwx_in;
 
@@ -190,7 +189,7 @@ void Paw_Cell::set_paw_k(
         {
             arg += atom_coord[iat][i] * kpt[i];
         }
-        arg *= twopi;
+        arg *= ModuleBase::TWO_PI;
         const std::complex<double> kphase = std::complex<double>(cos(arg), -sin(arg));
 
         struc_fact[iat].resize(npw);
@@ -326,17 +325,14 @@ std::vector<double> Paw_Cell::calc_ylm(const int lmax, const double * r)
     std::vector<std::complex<double>> phase;
     phase.resize(lmax+1);
 
-    const double pi = 3.141592653589793238462643383279502884197;
-    const double fourpi = 4.0 * pi;
-
     // set zero
     for(int i = 0; i < size_ylm; i++)
     {
         ylm[i] = 0;
     }
 
     // l = 0
-    ylm[0] = 1.0/std::sqrt(fourpi);
+    ylm[0] = 1.0 / std::sqrt(ModuleBase::FOUR_PI);
 
     if(rr>tol)
     {
@@ -361,8 +357,8 @@ std::vector<double> Paw_Cell::calc_ylm(const int lmax, const double * r)
         {
             const int l0 = l*l + l;
             double fact  = 1.0/double(l*(l+1));
-            const double ylmcst = std::sqrt(double(2*l+1)/fourpi);
-            
+            const double ylmcst = std::sqrt(double(2 * l + 1) / ModuleBase::FOUR_PI);
+
             // m = 0
             ylm[l0] = ylmcst*Paw_Cell::ass_leg_pol(l,0,ctheta);
 

source/module_cell/module_paw/paw_sphbes.cpp

@@ -1,11 +1,12 @@
-#include "paw_element.h"
-#include "module_base/tool_title.h"
+#include "module_base/constants.h"
 #include "module_base/tool_quit.h"
+#include "module_base/tool_title.h"
+#include "paw_element.h"
 
 double Paw_Element::get_ptilde(const int istate_in, const double q_in, const double omega)
 {
     // multiply by a factor 4pi / sqrt(omega)
-    const double factor = 4.0 * 3.141592653589793238462643383279502884197 / std::sqrt(omega);
+    const double factor = ModuleBase::FOUR_PI / std::sqrt(omega);
 
     return this->splint(qgrid, ptilde_q[istate_in], d2ptilde_q[istate_in], q_in) * factor;
 }
@@ -18,9 +19,6 @@ void Paw_Element::transform_ptilde()
     double dq = 0.01;
     int    nq = int( (std::sqrt(ecutwfc) / dq + 4) * cell_factor );
 
-    const double pi = 3.141592653589793238462643383279502884197;
-    const double twopi = 2.0 * pi;
-
     std::vector<double> integrand;
     integrand.resize(nr);
 
@@ -49,7 +47,7 @@ void Paw_Element::transform_ptilde()
         {
             for(int ir = 0; ir < nr; ir ++)
             {
-                integrand[ir] = twopi * ptilde_r[istate][ir] * std::pow(rr[ir],3);
+                integrand[ir] = ModuleBase::TWO_PI * ptilde_r[istate][ir] * std::pow(rr[ir], 3);
             }
             yp1 = this -> simpson_integration(integrand) / 3.0;
         }
@@ -60,7 +58,7 @@ void Paw_Element::transform_ptilde()
             double x = rr[ir] * double(nq - 1) * dq;
 
             this -> spherical_bessel_function(l,x,bes,besp,1);
-            integrand[ir] = twopi * besp * ptilde_r[istate][ir] * std::pow(rr[ir],3);
+            integrand[ir] = ModuleBase::TWO_PI * besp * ptilde_r[istate][ir] * std::pow(rr[ir], 3);
         }
 
         ypn = this -> simpson_integration(integrand);

source/module_elecstate/module_charge/charge_init.cpp

@@ -196,8 +196,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             {
                 std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
             }
-            const double pi = 3.141592653589790;
-            const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
+            const double fact = (3.0 / 5.0) * pow(3.0 * ModuleBase::PI * ModuleBase::PI, 2.0 / 3.0);
             for (int is = 0; is < PARAM.inp.nspin; ++is)
             {
                 for (int ir = 0; ir < this->rhopw->nrxx; ++ir)

source/module_hamilt_general/module_xc/xc_funct_exch_gga.cpp

@@ -96,8 +96,7 @@ double &sx, double &v1x, double &v2x)
 
 	// numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
     const double third = 1.0 / 3.0;
-    const double pi = 3.14159265358979323846;
-    const double c1 = 0.750 / pi;
+    const double c1 = 0.750 / ModuleBase::PI;
     const double c2 = 3.0936677262801360;
     const double c5 = 4.0 * third;
     // parameters of the functional
@@ -196,11 +195,10 @@ void XC_Functional::wcx(const double &rho,const double &grho, double &sx, double
   // exchange energy gradient part
   double dxunif, dfx, f1, f2, f3, dfx1, x1, x2, x3, dxds1, dxds2, dxds3;
   // numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
-  double third, c1, c2, c5, teneightyone;	// c6
-  const double pi = 3.14159265358979323846;
+  double third, c1, c2, c5, teneightyone; // c6
 
   third = 1.0/3.0;
-  c1 = 0.75/pi;
+  c1 = 0.75 / ModuleBase::PI;
   c2 = 3.093667726280136;
   c5 = 4.0 * third;
   teneightyone = 0.123456790123;

source/module_hamilt_general/module_xc/xc_funct_exch_lda.cpp

@@ -35,15 +35,14 @@ void XC_Functional::slater1(const double &rs, double &ex, double &vx)
 // Slater exchange with alpha=2/3 and Relativistic exchange
 void XC_Functional::slater_rxc(const double &rs, double &ex, double &vx)
 {
-	static const double pi = 3.14159265358979;
     const double trd = 1.0 / 3.0;
     //const double ftrd = 4.0 / 3.0;
     //const double tftm = pow(2.0, ftrd) - 2.0;
-    const double a0 = pow((4.0 / (9.0 * pi)), trd);
+    const double a0 = pow((4.0 / (9.0 * ModuleBase::PI)), trd);
     // X-alpha parameter:
     const double alp = 2 * trd;
 
-    double vxp = -3 * alp / (2 * pi * a0 * rs);
+    double vxp = -3 * alp / (2 * ModuleBase::PI * a0 * rs);
     double exp = 3 * vxp / 4;
     const double beta = 0.014 / rs;
     const double sb = sqrt(1 + beta * beta);
@@ -53,7 +52,7 @@ void XC_Functional::slater_rxc(const double &rs, double &ex, double &vx)
     exp = exp * (1.0 - 1.5 * x * x);
     vx = vxp;
     ex = exp;
-	return;	
+    return;
 }
 
 // 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,
         return;
     }
 
-	const double pi = 3.141592653589790;
-
-	const double trd = 1.0 / 3.0;
+    const double trd = 1.0 / 3.0;
 	const double ftrd = 4.0 / 3.0;
 	double tftm = pow(2.0, ftrd) - 2;
-	double a0 = pow((4 / (9 * pi)), trd);
+    double a0 = pow((4 / (9 * ModuleBase::PI)), trd);
 
-	double alp = 2 * trd;
+    double alp = 2 * trd;
 
 	double fz = (pow((1 + z), ftrd) + pow((1 - z), ftrd) - 2) / tftm;
 	double fzp = ftrd * (pow((1 + z), trd) - pow((1 - z), trd)) / tftm;
 
-    double rs = pow((3 / (4 * pi * rho)), trd);
-    double vxp = -3 * alp / (2 * pi * a0 * rs);
+    double rs = pow((3 / (4 * ModuleBase::PI * rho)), trd);
+    double vxp = -3 * alp / (2 * ModuleBase::PI * a0 * rs);
     double exp = 3 * vxp / 4;
 
     double beta = 0.014 / rs;

source/module_hamilt_general/module_xc/xc_funct_hcth.cpp

@@ -20,7 +20,6 @@ void XC_Functional::hcth(const double rho, const double grho, double &sx, double
     // real(kind=DP) :: rho, grho, sx, v1x, v2x
 
     // parameter :
-    const double pi = 3.14159265358979323846;
     double o3 = 1.00 / 3.00;
     double o34 = 4.00 / 3.00;
     double fr83 = 8.0 / 3.0;
@@ -33,7 +32,7 @@ void XC_Functional::hcth(const double rho, const double grho, double &sx, double
     dgaa_drho, dgab_drho, dgx_drho, dgaa_dgr, dgab_dgr, dgx_dgr;
 
     r3q2 = std::pow(2.0, (-o3));
-    r3pi = std::pow((3.0 / pi), o3);
+    r3pi = std::pow((3.0 / ModuleBase::PI), o3);
     //.....coefficients for pwf correlation......................................
     cg0[1] = 0.0310910;
     cg0[2] = 0.2137000;

source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp

@@ -176,19 +176,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
 
 	double *aux1 = new double[msh];
 
-	// for tests
-	/*
-	for(ir=0; ir<msh; ir++)
-	{
-		aux[ir] = r[ir] * zp_in * e2 / ucell.omega;
-	}
-	ModuleBase::Integral::Simpson_Integral(msh, aux, rab, vloc_1d[0] );
-	vloc_1d[0] *= 4*3.1415926;
-	std::cout << "  vloc_1d[0]=" <<  vloc_1d[0]/rho_basis->npw << std::endl;
-	std::cout << "  vloc_1d[0]=" <<  vloc_1d[0]/rho_basis->nxyz << std::endl;
-	*/
-
-	// (1)
+    // (1)
 	if(rho_basis->gg_uniq[0] < 1.0e-8)
 	{
 		double *aux = new double[msh];