diff --git a/source/module_basis/module_ao/ORB_nonlocal_lm.cpp b/source/module_basis/module_ao/ORB_nonlocal_lm.cpp
index ffb0fc795d6..9db01eb4a13 100644
--- a/source/module_basis/module_ao/ORB_nonlocal_lm.cpp
+++ b/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++)
     {
diff --git a/source/module_cell/module_neighbor/sltk_util.h b/source/module_cell/module_neighbor/sltk_util.h
index a9b14454226..0e783899c8b 100644
--- a/source/module_cell/module_neighbor/sltk_util.h
+++ b/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)
diff --git a/source/module_cell/module_paw/paw_cell.cpp b/source/module_cell/module_paw/paw_cell.cpp
index 7fd6337e15b..1d1e491c1d4 100644
--- a/source/module_cell/module_paw/paw_cell.cpp
+++ 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,9 +325,6 @@ 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++)
     {
@@ -336,7 +332,7 @@ std::vector<double> Paw_Cell::calc_ylm(const int lmax, const double * r)
     }
 
     // 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);
             
diff --git a/source/module_cell/module_paw/paw_sphbes.cpp b/source/module_cell/module_paw/paw_sphbes.cpp
index 71360ba34e4..85ec3040c73 100644
--- a/source/module_cell/module_paw/paw_sphbes.cpp
+++ b/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);
diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp
index 6e671f887aa..b42043cb968 100644
--- a/source/module_elecstate/module_charge/charge_init.cpp
+++ b/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)
diff --git a/source/module_hamilt_general/module_xc/xc_funct_exch_gga.cpp b/source/module_hamilt_general/module_xc/xc_funct_exch_gga.cpp
index 4045707085a..e1f1c0cf27e 100644
--- a/source/module_hamilt_general/module_xc/xc_funct_exch_gga.cpp
+++ b/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;
diff --git a/source/module_hamilt_general/module_xc/xc_funct_exch_lda.cpp b/source/module_hamilt_general/module_xc/xc_funct_exch_lda.cpp
index 14f32c5a47b..9d6cb6d955b 100644
--- a/source/module_hamilt_general/module_xc/xc_funct_exch_lda.cpp
+++ b/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;
diff --git a/source/module_hamilt_general/module_xc/xc_funct_hcth.cpp b/source/module_hamilt_general/module_xc/xc_funct_hcth.cpp
index 37c92418f82..d7c3333c99f 100644
--- a/source/module_hamilt_general/module_xc/xc_funct_hcth.cpp
+++ b/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;
diff --git a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp
index 486b866bf12..8ab3319ad36 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp
+++ b/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];
diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.cpp b/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.cpp
index 0cd0e1ab96b..8c26ca9f13d 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.cpp
@@ -1,12 +1,14 @@
 #include "module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h"
 
+#include "module_base/constants.h"
+#include "module_base/libm/libm.h"
 #include "module_base/math_polyint.h"
 #include "module_base/memory.h"
 #include "module_hamilt_pw/hamilt_pwdft/kernels/vnl_op.h"
 #include "vnl_tools.hpp"
-#include <vector>
+
 #include <iomanip>
-#include "module_base/libm/libm.h"
+#include <vector>
 namespace hamilt
 {
 template <typename FPTYPE>
@@ -567,9 +569,8 @@ struct cal_stress_drhoc_aux_op<FPTYPE, base_device::DEVICE_CPU> {
                     const int igl0,
                     const int ngg,
                     const double omega,
-                    int type) {
-    const double FOUR_PI = 4.0 * 3.14159265358979323846;
-    // printf("%d,%d,%lf\n",ngg,mesh,omega);
+                    int type)
+    {
 
 #ifdef _OPENMP
 #pragma omp parallel
@@ -601,13 +602,24 @@ struct cal_stress_drhoc_aux_op<FPTYPE, base_device::DEVICE_CPU> {
                 }
             }//ir
             Simpson_Integral<FPTYPE>(mesh, aux.data(), rab, rhocg1);
-            if(type ==0 ) { drhocg [igl] = FOUR_PI / omega * rhocg1;
-            } else if(type == 1) { drhocg [igl] = FOUR_PI * rhocg1 / omega;
-            } else if(type == 2) { drhocg [igl] = rhocg1;
-            } else if(type == 3) {
-                rhocg1 *= FOUR_PI / omega / 2.0 / gx_arr[igl];
+            if (type == 0)
+            {
+                drhocg[igl] = ModuleBase::FOUR_PI / omega * rhocg1;
+            }
+            else if (type == 1)
+            {
+                drhocg[igl] = ModuleBase::FOUR_PI * rhocg1 / omega;
+            }
+            else if (type == 2)
+            {
+                drhocg[igl] = rhocg1;
+            }
+            else if (type == 3)
+            {
+                rhocg1 *= ModuleBase::FOUR_PI / omega / 2.0 / gx_arr[igl];
                 FPTYPE g2a = (gx_arr[igl]*gx_arr[igl]) / 4.0;
-                rhocg1 += FOUR_PI / omega * gx_arr[ngg] * ModuleBase::libm::exp ( - g2a) * (g2a + 1) / pow(gx_arr[igl]*gx_arr[igl] , 2);
+                rhocg1 += ModuleBase::FOUR_PI / omega * gx_arr[ngg] * ModuleBase::libm::exp(-g2a) * (g2a + 1)
+                          / pow(gx_arr[igl] * gx_arr[igl], 2);
                 drhocg [igl] = rhocg1;
             }
         }
@@ -620,17 +632,20 @@ struct cal_stress_drhoc_aux_op<FPTYPE, base_device::DEVICE_CPU> {
 
 template <typename FPTYPE>
 struct cal_force_npw_op<FPTYPE, base_device::DEVICE_CPU> {
-    void operator()(const std::complex<FPTYPE> *psiv,
-                    const FPTYPE* gv_x, const FPTYPE* gv_y, const FPTYPE* gv_z,
+    void operator()(const std::complex<FPTYPE>* psiv,
+                    const FPTYPE* gv_x,
+                    const FPTYPE* gv_y,
+                    const FPTYPE* gv_z,
                     const FPTYPE* rhocgigg_vec,
                     FPTYPE* force,
-                    const FPTYPE pos_x, const FPTYPE pos_y, const FPTYPE pos_z,
+                    const FPTYPE pos_x,
+                    const FPTYPE pos_y,
+                    const FPTYPE pos_z,
                     const int npw,
-                    const FPTYPE omega, const FPTYPE tpiba) {
-        const double TWO_PI = 2.0 * 3.14159265358979323846;
-        // printf("%d,%d,%lf\n",ngg,mesh,omega);
+                    const FPTYPE omega,
+                    const FPTYPE tpiba)
+    {
 
-        // printf("iminininininin\n\n\n\n");
 #ifdef _OPENMP
 #pragma omp for nowait
 #endif
@@ -638,7 +653,7 @@ struct cal_force_npw_op<FPTYPE, base_device::DEVICE_CPU> {
         {
             const std::complex<FPTYPE> psiv_conj = conj(psiv[ig]);
 
-            const FPTYPE arg = TWO_PI * (gv_x[ig] * pos_x + gv_y[ig] * pos_y + gv_z[ig] * pos_z);
+            const FPTYPE arg = ModuleBase::TWO_PI * (gv_x[ig] * pos_x + gv_y[ig] * pos_y + gv_z[ig] * pos_z);
             FPTYPE sinp, cosp;
             ModuleBase::libm::sincos(arg, &sinp, &cosp);
             const std::complex<FPTYPE> expiarg = std::complex<FPTYPE>(sinp, cosp);
diff --git a/source/module_io/cal_dos.cpp b/source/module_io/cal_dos.cpp
index 6f18b12af96..f1383764a16 100644
--- a/source/module_io/cal_dos.cpp
+++ b/source/module_io/cal_dos.cpp
@@ -134,7 +134,7 @@ bool ModuleIO::calculate_dos(const int& is,
                 //----------------------------------------------------------
                 // EXPLAIN : if en
                 //----------------------------------------------------------
-                Gauss = exp(-de2 / b / b) / sqrt(3.1415926) / b;
+                Gauss = exp(-de2 / b / b) / sqrt(ModuleBase::PI) / b;
                 dos_smearing[j] += dos[i] * Gauss;
             }
         }
diff --git a/source/module_io/write_dos_lcao.cpp b/source/module_io/write_dos_lcao.cpp
index 14beb2830c5..4b4deafb66c 100644
--- a/source/module_io/write_dos_lcao.cpp
+++ b/source/module_io/write_dos_lcao.cpp
@@ -111,8 +111,7 @@ void ModuleIO::write_dos_lcao(const UnitCell& ucell,
     }
 
     double a = bcoeff;
-    double c = 2 * 3.141592653;
-    double b = sqrt(c) * a;
+    double b = sqrt(ModuleBase::TWO_PI) * a;
 
     std::complex<double>* waveg = new std::complex<double>[PARAM.globalv.nlocal];
 
@@ -428,8 +427,7 @@ void ModuleIO::write_dos_lcao(const UnitCell& ucell,
         }
 
         double a = bcoeff;
-        double c = 2 * 3.141592653;
-        double b = sqrt(c) * a;
+        double b = sqrt(ModuleBase::TWO_PI) * a;
 
         std::complex<double>* waveg = new std::complex<double>[PARAM.globalv.nlocal];
 
diff --git a/tests/integrate/204_NO_KP_NC/result.ref b/tests/integrate/204_NO_KP_NC/result.ref
index 7283435a661..933f82a8358 100644
--- a/tests/integrate/204_NO_KP_NC/result.ref
+++ b/tests/integrate/204_NO_KP_NC/result.ref
@@ -1,4 +1,4 @@
-etotref -6436.214785882292
-etotperatomref -3218.1073929411
+etotref -6436.214785807903
+etotperatomref -3218.1073929040
 totalforceref 0.000002
-totaltimeref 10.22
+totaltimeref 12.03
diff --git a/tests/integrate/218_NO_mixing_tau/result.ref b/tests/integrate/218_NO_mixing_tau/result.ref
index c6c8f33142d..b18a26208e0 100644
--- a/tests/integrate/218_NO_mixing_tau/result.ref
+++ b/tests/integrate/218_NO_mixing_tau/result.ref
@@ -1,3 +1,3 @@
-etotref -122.1473084278620
-etotperatomref -122.1473084279
-totaltimeref 2.60
+etotref -122.1473085261831
+etotperatomref -122.1473085262
+totaltimeref 2.88