diff --git a/src/CLASS2/angle_class2.cpp b/src/CLASS2/angle_class2.cpp
index ca69cca5b2..ac74f441f7 100644
--- a/src/CLASS2/angle_class2.cpp
+++ b/src/CLASS2/angle_class2.cpp
@@ -24,10 +24,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -318,7 +320,7 @@ void AngleClass2::coeff(int narg, char **arg)
     // convert theta0 from degrees to radians
 
     for (int i = ilo; i <= ihi; i++) {
-      theta0[i] = theta0_one/180.0 * PI;
+      theta0[i] = theta0_one/180.0 * MY_PI;
       k2[i] = k2_one;
       k3[i] = k3_one;
       k4[i] = k4_one;
diff --git a/src/CLASS2/dihedral_class2.cpp b/src/CLASS2/dihedral_class2.cpp
index d6e31696d7..f9c04126d8 100644
--- a/src/CLASS2/dihedral_class2.cpp
+++ b/src/CLASS2/dihedral_class2.cpp
@@ -26,20 +26,19 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 #define SMALL     0.0000001
 
 /* ---------------------------------------------------------------------- */
 
-DihedralClass2::DihedralClass2(LAMMPS *lmp) : Dihedral(lmp)
-{
-  PI = 4.0*atan(1.0);
-}
+DihedralClass2::DihedralClass2(LAMMPS *lmp) : Dihedral(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
@@ -697,8 +696,8 @@ void DihedralClass2::coeff(int narg, char **arg)
       at_f1_2[i] = f1_2_one;
       at_f2_2[i] = f2_2_one;
       at_f3_2[i] = f3_2_one;
-      at_theta0_1[i] = theta0_1_one/180.0 * PI;
-      at_theta0_2[i] = theta0_2_one/180.0 * PI;
+      at_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
+      at_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
       setflag_at[i] = 1;
       count++;
     }
@@ -714,8 +713,8 @@ void DihedralClass2::coeff(int narg, char **arg)
     
     for (int i = ilo; i <= ihi; i++) {
       aat_k[i] = k_one;
-      aat_theta0_1[i] = theta0_1_one/180.0 * PI;
-      aat_theta0_2[i] = theta0_2_one/180.0 * PI;
+      aat_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
+      aat_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
       setflag_aat[i] = 1;
       count++;
     }
@@ -749,11 +748,11 @@ void DihedralClass2::coeff(int narg, char **arg)
 
     for (int i = ilo; i <= ihi; i++) {
       k1[i] = k1_one;
-      phi1[i] = phi1_one/180.0 * PI;
+      phi1[i] = phi1_one/180.0 * MY_PI;
       k2[i] = k2_one;
-      phi2[i] = phi2_one/180.0 * PI;
+      phi2[i] = phi2_one/180.0 * MY_PI;
       k3[i] = k3_one;
-      phi3[i] = phi3_one/180.0 * PI;
+      phi3[i] = phi3_one/180.0 * MY_PI;
       setflag_d[i] = 1;
       count++;
     }
diff --git a/src/CLASS2/dihedral_class2.h b/src/CLASS2/dihedral_class2.h
index c99258a627..b7a519a27d 100644
--- a/src/CLASS2/dihedral_class2.h
+++ b/src/CLASS2/dihedral_class2.h
@@ -46,7 +46,6 @@ class DihedralClass2 : public Dihedral {
   double *bb13t_k,*bb13t_r10,*bb13t_r30;
   int *setflag_d,*setflag_mbt,*setflag_ebt;
   int *setflag_at,*setflag_aat,*setflag_bb13t;
-  double PI;
 
   void allocate();
 };
diff --git a/src/CLASS2/improper_class2.cpp b/src/CLASS2/improper_class2.cpp
index 7329870aec..c1331c607e 100644
--- a/src/CLASS2/improper_class2.cpp
+++ b/src/CLASS2/improper_class2.cpp
@@ -26,19 +26,18 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
 /* ---------------------------------------------------------------------- */
 
-ImproperClass2::ImproperClass2(LAMMPS *lmp) : Improper(lmp)
-{
-  PI = 4.0*atan(1.0);
-}
+ImproperClass2::ImproperClass2(LAMMPS *lmp) : Improper(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
@@ -550,9 +549,9 @@ void ImproperClass2::coeff(int narg, char **arg)
       aa_k1[i] = k1_one;
       aa_k2[i] = k2_one;
       aa_k3[i] = k3_one;
-      aa_theta0_1[i] = theta0_1_one/180.0 * PI;
-      aa_theta0_2[i] = theta0_2_one/180.0 * PI;
-      aa_theta0_3[i] = theta0_3_one/180.0 * PI;
+      aa_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
+      aa_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
+      aa_theta0_3[i] = theta0_3_one/180.0 * MY_PI;
       setflag_aa[i] = 1;
       count++;
     }
@@ -567,7 +566,7 @@ void ImproperClass2::coeff(int narg, char **arg)
 
     for (int i = ilo; i <= ihi; i++) {
       k0[i] = k0_one;
-      chi0[i] = chi0_one/180.0 * PI;
+      chi0[i] = chi0_one/180.0 * MY_PI;
       setflag_i[i] = 1;
       count++;
     }
diff --git a/src/CLASS2/improper_class2.h b/src/CLASS2/improper_class2.h
index 56c7c9d14a..9a86e9cd86 100644
--- a/src/CLASS2/improper_class2.h
+++ b/src/CLASS2/improper_class2.h
@@ -38,7 +38,6 @@ class ImproperClass2 : public Improper {
   double *k0,*chi0;
   double *aa_k1,*aa_k2,*aa_k3,*aa_theta0_1,*aa_theta0_2,*aa_theta0_3;
   int *setflag_i,*setflag_aa;
-  double PI;
 
   void allocate();
   void angleangle(int, int);
diff --git a/src/CLASS2/pair_lj_class2.cpp b/src/CLASS2/pair_lj_class2.cpp
index aba3715cc1..ee5d7d3e32 100644
--- a/src/CLASS2/pair_lj_class2.cpp
+++ b/src/CLASS2/pair_lj_class2.cpp
@@ -19,10 +19,12 @@
 #include "comm.h"
 #include "force.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -254,14 +256,13 @@ double PairLJClass2::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
     double sig6 = sig3*sig3;
     double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
     double rc6 = rc3*rc3;
-    etail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] *
+    etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
       sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
-    ptail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig3 - 2.0*rc3) / rc6;
   } 
 
diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp
index b6dc7d8018..9cd784ff58 100644
--- a/src/CLASS2/pair_lj_class2_coul_cut.cpp
+++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp
@@ -20,10 +20,12 @@
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -307,14 +309,13 @@ double PairLJClass2CoulCut::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
     double sig6 = sig3*sig3;
     double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
     double rc6 = rc3*rc3;
-    etail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] *
+    etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
       sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
-    ptail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig3 - 2.0*rc3) / rc6;
   } 
 
diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp
index b9423c4728..34ff47d804 100644
--- a/src/CLASS2/pair_lj_class2_coul_long.cpp
+++ b/src/CLASS2/pair_lj_class2_coul_long.cpp
@@ -22,10 +22,12 @@
 #include "kspace.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define EWALD_F   1.12837917
 #define EWALD_P   0.3275911
@@ -321,14 +323,13 @@ double PairLJClass2CoulLong::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
     double sig6 = sig3*sig3;
     double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
     double rc6 = rc3*rc3;
-    etail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] *
+    etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
       sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
-    ptail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig3 - 2.0*rc3) / rc6;
   } 
 
diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp
index 5eb072b93e..f3c9a0dc06 100644
--- a/src/COLLOID/pair_lubricate.cpp
+++ b/src/COLLOID/pair_lubricate.cpp
@@ -28,11 +28,13 @@
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "update.h"
-#include "memory.h"
 #include "random_mars.h"
+#include "math_const.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -73,7 +75,6 @@ void PairLubricate::compute(int eflag, int vflag)
   double P_dot_wrel_1,P_dot_wrel_2,P_dot_wrel_3;
   double a_squeeze,a_shear,a_pump,a_twist;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  double PI = 4.0*atan(1.0);
 
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
@@ -190,16 +191,16 @@ void PairLubricate::compute(int eflag, int vflag)
         h_sep = r - 2.0*radi;
 
 	if (flag1)
-	  a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep);
+	  a_squeeze = (3.0*MY_PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep);
 	if (flag2) 
-	  a_shear = (PI*mu*2.*radi/2.0) *
+	  a_shear = (MY_PI*mu*2.*radi/2.0) *
 	    log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0;
 	if (flag3) 
-	  a_pump = (PI*mu*pow(2.0*radi,4)/8.0) *
+	  a_pump = (MY_PI*mu*pow(2.0*radi,4)/8.0) *
 	    ((3.0/20.0) * log(2.0*radi/2.0/h_sep) + 
 	     (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep));
 	if (flag4)
-	  a_twist = (PI*mu*pow(2.0*radi,4)/4.0) *
+	  a_twist = (MY_PI*mu*pow(2.0*radi,4)/4.0) *
 	    (h_sep/2.0/radi) * log(2.0/(2.0*h_sep));
 
         if (h_sep >= cut_inner[itype][jtype]) {
@@ -231,7 +232,7 @@ void PairLubricate::compute(int eflag, int vflag)
 	  torque[i][2] += vxmu2f * tz;
 
         } else {
-	  a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * 
+	  a_squeeze = (3.0*MY_PI*mu*2.0*radi/2.0) * 
 	    (2.0*radi/4.0/cut_inner[itype][jtype]);
 	  fpair = -a_squeeze*vnnr;
 	  fpair *= vxmu2f;
diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp
index 044d018c28..0fa91a2518 100644
--- a/src/GPU/pppm_gpu.cpp
+++ b/src/GPU/pppm_gpu.cpp
@@ -32,11 +32,13 @@
 #include "domain.h"
 #include "fft3d_wrap.h"
 #include "remap_wrap.h"
+#include "gpu_extra.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
-#include "gpu_extra.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define MAXORDER 7
 #define OFFSET 16384
@@ -189,7 +191,7 @@ void PPPMGPU::compute(int eflag, int vflag)
    
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/1.772453851 +
-      0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
+      MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e*scale;
   }
 
diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 4dcf9fc3bc..03a1c3a62a 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -27,10 +27,12 @@
 #include "region_block.h"
 #include "region_cylinder.h"
 #include "random_park.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define EPSILON 0.001
 
@@ -54,8 +56,6 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
 
   if (seed <= 0) error->all(FLERR,"Illegal fix pour command");
 
-  PI = 4.0*atan(1.0);
-
   // option defaults
 
   int iregion = -1;
@@ -214,11 +214,11 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
       double dy = yhi - ylo;
       if (dy < 1.0) dy = 1.0;
       volume = (xhi-xlo) * dy * (zhi-zlo);
-    } else volume = PI*rc*rc * (zhi-zlo);
-    volume_one = 4.0/3.0 * PI * radius_hi*radius_hi*radius_hi;
+    } else volume = MY_PI*rc*rc * (zhi-zlo);
+    volume_one = 4.0/3.0 * MY_PI * radius_hi*radius_hi*radius_hi;
   } else {
     volume = (xhi-xlo) * (yhi-ylo);
-    volume_one = PI * radius_hi*radius_hi;
+    volume_one = MY_PI * radius_hi*radius_hi;
   }
 
   nper = static_cast<int> (volfrac*volume/volume_one);
@@ -472,7 +472,7 @@ void FixPour::pre_exchange()
       m = atom->nlocal - 1;
       atom->type[m] = ntype;
       atom->radius[m] = radtmp;
-      atom->rmass[m] = 4.0*PI/3.0 * radtmp*radtmp*radtmp * denstmp;
+      atom->rmass[m] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp;
       atom->mask[m] = 1 | groupbit;
       atom->v[m][0] = vxtmp;
       atom->v[m][1] = vytmp;
diff --git a/src/GRANULAR/fix_pour.h b/src/GRANULAR/fix_pour.h
index ec9590b276..1880597139 100644
--- a/src/GRANULAR/fix_pour.h
+++ b/src/GRANULAR/fix_pour.h
@@ -56,7 +56,6 @@ class FixPour : public Fix {
 
   int me,nprocs;
   int *recvcounts,*displs;
-  double PI;
   int nfreq,nfirst,ninserted,nper;
   double lo_current,hi_current;
   class FixShearHistory *fix_history;
diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp
index 0345274cb0..020425d110 100644
--- a/src/GRANULAR/fix_wall_gran.cpp
+++ b/src/GRANULAR/fix_wall_gran.cpp
@@ -26,10 +26,12 @@
 #include "pair.h"
 #include "modify.h"
 #include "respa.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{XPLANE,YPLANE,ZPLANE,ZCYLINDER};    // XYZ PLANE need to be 0,1,2
 enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY};
@@ -159,10 +161,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
 
   // setup oscillations
 
-  if (wiggle) {
-    double PI = 4.0 * atan(1.0);
-    omega = 2.0*PI / period;
-  }
+  if (wiggle) omega = 2.0*MY_PI / period;
 
   // perform initial allocation of atom-based arrays
   // register with Atom class
diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index 5b6d46750a..edd030b55b 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -26,10 +26,12 @@
 #include "force.h"
 #include "pair.h"
 #include "domain.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.00001
 
@@ -40,7 +42,6 @@ Ewald::Ewald(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
   if (narg != 1) error->all(FLERR,"Illegal kspace_style ewald command");
 
   precision = atof(arg[0]);
-  PI = 4.0*atan(1.0);
 
   kmax = 0;
   kxvecs = kyvecs = kzvecs = NULL;
@@ -165,17 +166,17 @@ void Ewald::setup()
   double zprd_slab = zprd*slab_volfactor;
   volume = xprd * yprd * zprd_slab;
 
-  unitk[0] = 2.0*PI/xprd;
-  unitk[1] = 2.0*PI/yprd;
-  unitk[2] = 2.0*PI/zprd_slab;
+  unitk[0] = 2.0*MY_PI/xprd;
+  unitk[1] = 2.0*MY_PI/yprd;
+  unitk[2] = 2.0*MY_PI/zprd_slab;
 
   // determine kmax
   // function of current box size, precision, G_ewald (short-range cutoff)
 
-  int nkxmx = static_cast<int> ((g_ewald*xprd/PI) * sqrt(-log(precision)));
-  int nkymx = static_cast<int> ((g_ewald*yprd/PI) * sqrt(-log(precision)));
+  int nkxmx = static_cast<int> ((g_ewald*xprd/MY_PI) * sqrt(-log(precision)));
+  int nkymx = static_cast<int> ((g_ewald*yprd/MY_PI) * sqrt(-log(precision)));
   int nkzmx = 
-    static_cast<int> ((g_ewald*zprd_slab/PI) * sqrt(-log(precision)));
+    static_cast<int> ((g_ewald*zprd_slab/MY_PI) * sqrt(-log(precision)));
 
   int kmax_old = kmax;
   kmax = MAX(nkxmx,nkymx);
@@ -281,9 +282,8 @@ void Ewald::compute(int eflag, int vflag)
     for (k = 0; k < kcount; k++)
       energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + 
 			 sfacim_all[k]*sfacim_all[k]);
-    PI = 4.0*atan(1.0);
     energy -= g_ewald*qsqsum/1.772453851 + 
-      0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
+      MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e*scale;
   }
 
@@ -495,7 +495,7 @@ void Ewald::coeffs()
   double unitky = unitk[1];
   double unitkz = unitk[2];
   double g_ewald_sq_inv = 1.0 / (g_ewald*g_ewald);
-  double preu = 4.0*PI/volume;
+  double preu = 4.0*MY_PI/volume;
 
   kcount = 0;
 
@@ -817,13 +817,13 @@ void Ewald::slabcorr(int eflag)
 
   // compute corrections
   
-  double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume;
+  double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
   
   if (eflag) energy += qqrd2e*scale * e_slabcorr;
 
   // add on force corrections
 
-  double ffact = -4.0*PI*dipole_all/volume; 
+  double ffact = -4.0*MY_PI*dipole_all/volume; 
   double **f = atom->f;
 
   for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;
diff --git a/src/KSPACE/ewald.h b/src/KSPACE/ewald.h
index 250b524de0..95482cfbfe 100644
--- a/src/KSPACE/ewald.h
+++ b/src/KSPACE/ewald.h
@@ -34,7 +34,6 @@ class Ewald : public KSpace {
   double memory_usage();
 
  private:
-  double PI;
   double precision;
   int kcount,kmax,kmax3d,kmax_created;
   double qqrd2e;
diff --git a/src/KSPACE/pair_born_coul_long.cpp b/src/KSPACE/pair_born_coul_long.cpp
index dca8532d35..2e607cd74d 100644
--- a/src/KSPACE/pair_born_coul_long.cpp
+++ b/src/KSPACE/pair_born_coul_long.cpp
@@ -26,10 +26,12 @@
 #include "kspace.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define EWALD_F   1.12837917
 #define EWALD_P   0.3275911
@@ -313,7 +315,6 @@ double PairBornCoulLong::init_one(int i, int j)
      } 
      MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
 
-     double PI = 4.0*atan(1.0);
      double rho1 = rho[i][j];
      double rho2 = rho1*rho1;
      double rho3 = rho2*rho1;
@@ -321,11 +322,11 @@ double PairBornCoulLong::init_one(int i, int j)
      double rc2 = rc*rc;
      double rc3 = rc2*rc;
      double rc5 = rc3*rc2;
-     etail_ij = 2.0*PI*all[0]*all[1] * 
+     etail_ij = 2.0*MY_PI*all[0]*all[1] * 
        (a[i][j]*exp((sigma[i][j]-rc)/rho1)*rho1* 
 	(rc2 + 2.0*rho1*rc + 2.0*rho2) - 
 	c[i][j]/(3.0*rc3) + d[i][j]/(5.0*rc5));
-     ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1] * 
+     ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1] * 
        (-a[i][j]*exp((sigma[i][j]-rc)/rho1) * 
 	(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 
 	2.0*c[i][j]/rc3 - 8.0*d[i][j]/(5.0*rc5)); 
diff --git a/src/KSPACE/pair_buck_coul_long.cpp b/src/KSPACE/pair_buck_coul_long.cpp
index 813b3b025a..a68029f367 100644
--- a/src/KSPACE/pair_buck_coul_long.cpp
+++ b/src/KSPACE/pair_buck_coul_long.cpp
@@ -22,10 +22,12 @@
 #include "kspace.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define EWALD_F   1.12837917
 #define EWALD_P   0.3275911
@@ -293,17 +295,16 @@ double PairBuckCoulLong::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
 
-    double PI = 4.0*atan(1.0);
     double rho1 = rho[i][j];
     double rho2 = rho1*rho1;
     double rho3 = rho2*rho1;
     double rc = cut_lj[i][j];
     double rc2 = rc*rc;
     double rc3 = rc2*rc;
-    etail_ij = 2.0*PI*all[0]*all[1]*
+    etail_ij = 2.0*MY_PI*all[0]*all[1]*
       (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) - 
        c[i][j]/(3.0*rc3));
-    ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1]*
+    ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1]*
       (-a[i][j]*exp(-rc/rho1)*
        (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3);
   }
diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp
index ba9e390baa..ff8fb1fd22 100644
--- a/src/KSPACE/pair_lj_cut_coul_long.cpp
+++ b/src/KSPACE/pair_lj_cut_coul_long.cpp
@@ -30,10 +30,12 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define EWALD_F   1.12837917
 #define EWALD_P   0.3275911
@@ -770,15 +772,14 @@ double PairLJCutCoulLong::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig2 = sigma[i][j]*sigma[i][j];
     double sig6 = sig2*sig2*sig2;
     double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
     double rc6 = rc3*rc3;
     double rc9 = rc3*rc6;
-    etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); 
-    ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); 
   } 
 
diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp
index e8c623b47b..3f0cb6e253 100644
--- a/src/MANYBODY/pair_airebo.cpp
+++ b/src/MANYBODY/pair_airebo.cpp
@@ -29,10 +29,12 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define MAXLINE 1024
 #define TOL 1.0e-9
@@ -52,8 +54,6 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp)
   maxpage = 0;
   pages = NULL;
   nC = nH = NULL;
-
-  PI = 4.0*atan(1.0);
 }
 
 /* ----------------------------------------------------------------------
@@ -1218,8 +1218,8 @@ double PairAIREBO::Sp(double Xij, double Xmin, double Xmax, double &dX)
     dX = 0.0;
   } 
   else {
-    cutoff = 0.5 * (1.0+cos(PI*t));
-    dX = (-0.5*PI*sin(PI*t)) / (Xmax-Xmin);
+    cutoff = 0.5 * (1.0+cos(MY_PI*t));
+    dX = (-MY_PI2*sin(MY_PI*t)) / (Xmax-Xmin);
   }
   return cutoff;
 }
diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h
index b51c00bf2b..166e1ed4ed 100644
--- a/src/MANYBODY/pair_airebo.h
+++ b/src/MANYBODY/pair_airebo.h
@@ -46,7 +46,6 @@ class PairAIREBO : public Pair {
   int npage;                       // current page in page list
   int *map;                        // 0 (C), 1 (H), or -1 (NULL) for each type
 
-  double PI;
   double cutlj;                    // user-specified LJ cutoff
   double cutljrebosq;              // cut for when to compute
                                    // REBO neighs of ghost atoms
diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp
index 5c47b8dcf1..2311ea8d1e 100644
--- a/src/MANYBODY/pair_comb.cpp
+++ b/src/MANYBODY/pair_comb.cpp
@@ -29,12 +29,14 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
-#include "memory.h"
-#include "error.h"
 #include "group.h"
 #include "update.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define MAXLINE 1024
 #define DELTA 4
@@ -47,11 +49,6 @@ PairComb::PairComb(LAMMPS *lmp) : Pair(lmp)
   single_enable = 0;
   one_coeff = 1;
 
-  PI = 4.0*atan(1.0);
-  PI2 = 2.0*atan(1.0);
-  PI4 = atan(1.0);
-  PIsq = sqrt(PI);
-
   nmax = 0;
   NCo = NULL;
   bbij = NULL;
@@ -930,7 +927,7 @@ double PairComb::elp(Param *param, double rsqij, double rsqik,
     double rij,rik,costheta,lp1,lp3,lp6;
     double rmu,rmu2,comtt,fck;
     double pplp1 = param->plp1, pplp3 = param->plp3, pplp6 = param->plp6;
-    double c123 = cos(param->a123*PI/180.0);
+    double c123 = cos(param->a123*MY_PI/180.0);
     
     // cos(theta) of the i-j-k
     // cutoff function of rik
@@ -984,7 +981,7 @@ void PairComb::flp(Param *param, double rsqij, double rsqik,
     double pplp1 = param->plp1;
     double pplp3 = param->plp3;
     double pplp6 = param->plp6;
-    double c123 = cos(param->a123*PI/180.0);
+    double c123 = cos(param->a123*MY_PI/180.0);
     
     // fck_d = derivative of cutoff function
     
@@ -1083,7 +1080,7 @@ double PairComb::comb_fc(double r, Param *param)
   
   if (r < comb_R-comb_D) return 1.0;
   if (r > comb_R+comb_D) return 0.0;
-  return 0.5*(1.0 + cos(PI*(r - comb_R)/comb_D));
+  return 0.5*(1.0 + cos(MY_PI*(r - comb_R)/comb_D));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1095,7 +1092,7 @@ double PairComb::comb_fc_d(double r, Param *param)
   
   if (r < comb_R-comb_D) return 0.0;
   if (r > comb_R+comb_D) return 0.0;
-  return -(PI2/comb_D) * sin(PI*(r - comb_R)/comb_D);
+  return -(MY_PI2/comb_D) * sin(MY_PI*(r - comb_R)/comb_D);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1107,7 +1104,7 @@ double PairComb::comb_fc2(double r)
   
   if (r < comb_R) return 0.0;
   if (r > comb_D) return 1.0;
-  return 0.5*(1.0 + cos(PI*(r - comb_R)/(comb_D-comb_R)));
+  return 0.5*(1.0 + cos(MY_PI*(r - comb_R)/(comb_D-comb_R)));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1119,7 +1116,7 @@ double PairComb::comb_fc2_d(double r)
   
   if (r < comb_R) return 0.0;
   if (r > comb_D) return 0.0;
-  return -(PI2/(comb_D-comb_R)) * sin(PI*(r - comb_R)/(comb_D-comb_R));
+  return -(MY_PI2/(comb_D-comb_R)) * sin(MY_PI*(r - comb_R)/(comb_D-comb_R));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1131,7 +1128,7 @@ double PairComb::comb_fc3(double r)
   
   if (r < comb_R) return 1.0;
   if (r > comb_D) return 0.0;
-  return 0.5*(1.0 + cos(PI*(r - comb_R)/(comb_D-comb_R)));
+  return 0.5*(1.0 + cos(MY_PI*(r - comb_R)/(comb_D-comb_R)));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1143,7 +1140,7 @@ double PairComb::comb_fc3_d(double r)
   
   if (r < comb_R) return 0.0;
   if (r > comb_D) return 0.0;
-  return -(PI2/(comb_D-comb_R)) * sin(PI*(r - comb_R)/(comb_D-comb_R));
+  return -(MY_PI2/(comb_D-comb_R)) * sin(MY_PI*(r - comb_R)/(comb_D-comb_R));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1541,10 +1538,10 @@ void PairComb::potal_calc(double &calc1, double &calc2, double &calc3)
   alf = 0.20;
   esucon = force->qqr2e;
 
-  calc2 = (erfc(rcoul*alf)/rcoul/rcoul+2.0*alf/PIsq*
+  calc2 = (erfc(rcoul*alf)/rcoul/rcoul+2.0*alf/MY_PIS*
 	   exp(-alf*alf*rcoul*rcoul)/rcoul)*esucon/rcoul;
   calc3 = (erfc(rcoul*alf)/rcoul)*esucon;
-  calc1 = -(alf/PIsq*esucon+calc3*0.5);
+  calc1 = -(alf/MY_PIS*esucon+calc3*0.5);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1590,7 +1587,7 @@ void PairComb::direct(int inty, int mr1, int mr2, int mr3, double rsq,
  r = sqrt(rsq);
  r3 = r * rsq;
  alf = 0.20;
- alfdpi = 2.0*alf/PIsq;
+ alfdpi = 2.0*alf/MY_PIS;
  esucon = force->qqr2e;
  pot_tmp = 0.0;
  pot_d = 0.0;
diff --git a/src/MANYBODY/pair_comb.h b/src/MANYBODY/pair_comb.h
index c8782b4bbf..bda2805d56 100644
--- a/src/MANYBODY/pair_comb.h
+++ b/src/MANYBODY/pair_comb.h
@@ -58,7 +58,6 @@ class PairComb : public Pair {
     int powermint;
   };
   
-  double PI,PI2,PI4,PIsq;
   double cutmax;                // max cutoff for all elements
   int nelements;                // # of unique elements
   char **elements;              // names of unique elements
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 9016a26b28..9c8654100d 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -30,10 +30,12 @@
 #include "random_park.h"
 #include "force.h"
 #include "pair.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -71,10 +73,9 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   // compute beta, lambda, sigma, and the zz factor
   
   beta = 1.0/(force->boltz*reservoir_temperature);
-  double PI = 4.0*atan(1.0);
   double gas_mass = atom->mass[ntype];
   double lambda = sqrt(force->hplanck*force->hplanck/
-		       (2.0*PI*gas_mass*force->mvv2e*
+		       (2.0*MY_PI*gas_mass*force->mvv2e*
 			force->boltz*reservoir_temperature));
   sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e);
   zz = exp(beta*chemical_potential)/(pow(lambda,3));
diff --git a/src/MOLECULE/angle_charmm.cpp b/src/MOLECULE/angle_charmm.cpp
index 5dd81580cc..0361209e66 100644
--- a/src/MOLECULE/angle_charmm.cpp
+++ b/src/MOLECULE/angle_charmm.cpp
@@ -23,10 +23,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -209,7 +211,7 @@ void AngleCharmm::coeff(int narg, char **arg)
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = k_one;
-    theta0[i] = theta0_one/180.0 * PI;
+    theta0[i] = theta0_one/180.0 * MY_PI;
     k_ub[i] = k_ub_one;
     r_ub[i] = r_ub_one;
     setflag[i] = 1;
diff --git a/src/MOLECULE/angle_cosine.cpp b/src/MOLECULE/angle_cosine.cpp
index 6c6fa4fc87..dc82378083 100644
--- a/src/MOLECULE/angle_cosine.cpp
+++ b/src/MOLECULE/angle_cosine.cpp
@@ -19,10 +19,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -174,7 +176,7 @@ void AngleCosine::coeff(int narg, char **arg)
 
 double AngleCosine::equilibrium_angle(int i)
 {
-  return PI;
+  return MY_PI;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MOLECULE/angle_cosine_periodic.cpp b/src/MOLECULE/angle_cosine_periodic.cpp
index d401c1fffc..8aaf54cb9f 100644
--- a/src/MOLECULE/angle_cosine_periodic.cpp
+++ b/src/MOLECULE/angle_cosine_periodic.cpp
@@ -23,10 +23,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -226,7 +228,7 @@ void AngleCosinePeriodic::coeff(int narg, char **arg)
 
 double AngleCosinePeriodic::equilibrium_angle(int i)
 {
-  return PI;
+  return MY_PI;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MOLECULE/angle_cosine_squared.cpp b/src/MOLECULE/angle_cosine_squared.cpp
index 9262700679..95c1a3b922 100644
--- a/src/MOLECULE/angle_cosine_squared.cpp
+++ b/src/MOLECULE/angle_cosine_squared.cpp
@@ -23,10 +23,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -178,7 +180,7 @@ void AngleCosineSquared::coeff(int narg, char **arg)
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = k_one;
-    theta0[i] = theta0_one/180.0 * PI;
+    theta0[i] = theta0_one/180.0 * MY_PI;
     setflag[i] = 1;
     count++;
   }
diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp
index e3def97499..8a462abaed 100644
--- a/src/MOLECULE/angle_harmonic.cpp
+++ b/src/MOLECULE/angle_harmonic.cpp
@@ -19,10 +19,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -178,7 +180,7 @@ void AngleHarmonic::coeff(int narg, char **arg)
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = k_one;
-    theta0[i] = theta0_one/180.0 * PI;
+    theta0[i] = theta0_one/180.0 * MY_PI;
     setflag[i] = 1;
     count++;
   }
diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp
index 4395c20413..a8e6254500 100644
--- a/src/MOLECULE/angle_table.cpp
+++ b/src/MOLECULE/angle_table.cpp
@@ -24,10 +24,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{LINEAR,SPLINE};
 
@@ -238,8 +240,8 @@ void AngleTable::coeff(int narg, char **arg)
   // convert theta from degrees to radians
 
   for (int i = 0; i < tb->ninput; i++){
-    tb->afile[i] *= PI/180.0;
-    tb->ffile[i] *= 180.0/PI; 
+    tb->afile[i] *= MY_PI/180.0;
+    tb->ffile[i] *= 180.0/MY_PI; 
   }
 
   // spline read-in and compute a,e,f vectors within table
@@ -442,7 +444,7 @@ void AngleTable::compute_table(Table *tb)
   // delta = table spacing in angle for N-1 bins
 
   int tlm1 = tablength-1;
-  tb->delta = PI/ tlm1;
+  tb->delta = MY_PI / tlm1;
   tb->invdelta = 1.0/tb->delta;
   tb->deltasq6 = tb->delta*tb->delta / 6.0;
   
@@ -501,8 +503,8 @@ void AngleTable::param_extract(Table *tb, char *line)
       tb->fplo = atof(word);
       word = strtok(NULL," \t\n\r\f");
       tb->fphi = atof(word);
-      tb->fplo *= (180.0/PI)*(180.0/PI);
-      tb->fphi *= (180.0/PI)*(180.0/PI);
+      tb->fplo *= (180.0/MY_PI)*(180.0/MY_PI);
+      tb->fphi *= (180.0/MY_PI)*(180.0/MY_PI);
     } else if (strcmp(word,"EQ") == 0) {
       word = strtok(NULL," \t\n\r\f");
       tb->theta0 = atof(word);
diff --git a/src/MOLECULE/dihedral_charmm.cpp b/src/MOLECULE/dihedral_charmm.cpp
index a404311150..c762a56ad5 100644
--- a/src/MOLECULE/dihedral_charmm.cpp
+++ b/src/MOLECULE/dihedral_charmm.cpp
@@ -27,10 +27,12 @@
 #include "force.h"
 #include "pair.h"
 #include "update.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 
@@ -344,14 +346,12 @@ void DihedralCharmm::coeff(int narg, char **arg)
   if (weight_one < 0.0 || weight_one > 1.0) 
     error->all(FLERR,"Incorrect weight arg for dihedral coefficients");
 
-  double PI = 4.0*atan(1.0);
-                       
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = k_one;
     shift[i] = shift_one;
-    cos_shift[i] = cos(PI*shift_one/180.0);
-    sin_shift[i] = sin(PI*shift_one/180.0);
+    cos_shift[i] = cos(MY_PI*shift_one/180.0);
+    sin_shift[i] = sin(MY_PI*shift_one/180.0);
     multiplicity[i] = multiplicity_one;
     weight[i] = weight_one;
     setflag[i] = 1;
@@ -420,10 +420,9 @@ void DihedralCharmm::read_restart(FILE *fp)
   MPI_Bcast(&shift[1],atom->ndihedraltypes,MPI_INT,0,world);
   MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
 
-  double PI = 4.0*atan(1.0);
   for (int i = 1; i <= atom->ndihedraltypes; i++) {
     setflag[i] = 1;
-    cos_shift[i] = cos(PI*shift[i]/180.0);
-    sin_shift[i] = sin(PI*shift[i]/180.0);
+    cos_shift[i] = cos(MY_PI*shift[i]/180.0);
+    sin_shift[i] = sin(MY_PI*shift[i]/180.0);
   }
 }
diff --git a/src/MOLECULE/dihedral_helix.cpp b/src/MOLECULE/dihedral_helix.cpp
index 7a50eb4fc0..ee93d2a08d 100644
--- a/src/MOLECULE/dihedral_helix.cpp
+++ b/src/MOLECULE/dihedral_helix.cpp
@@ -27,10 +27,12 @@
 #include "comm.h"
 #include "force.h"
 #include "update.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 #define SMALL     0.001
@@ -192,9 +194,9 @@ void DihedralHelix::compute(int eflag, int vflag)
     siinv = 1.0/si;
 
     p = aphi[type]*(1.0 - c) + bphi[type]*(1.0 + cos(3.0*phi)) +
-      cphi[type]*(1.0 + cos(phi + 0.25*PI));
+      cphi[type]*(1.0 + cos(phi + MY_PI4));
     pd = -aphi[type] + 3.0*bphi[type]*sin(3.0*phi)*siinv +
-      cphi[type]*sin(phi + 0.25*PI)*siinv;
+      cphi[type]*sin(phi + MY_PI4)*siinv;
 
     if (eflag) edihedral = p;
 
diff --git a/src/MOLECULE/improper_harmonic.cpp b/src/MOLECULE/improper_harmonic.cpp
index 62d5cd54b6..87d94079bd 100644
--- a/src/MOLECULE/improper_harmonic.cpp
+++ b/src/MOLECULE/improper_harmonic.cpp
@@ -22,10 +22,12 @@
 #include "domain.h"
 #include "force.h"
 #include "update.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 #define SMALL     0.001
@@ -248,7 +250,7 @@ void ImproperHarmonic::coeff(int narg, char **arg)
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = k_one;
-    chi[i] = chi_one/180.0 * PI;
+    chi[i] = chi_one/180.0 * MY_PI;
     setflag[i] = 1;
     count++;
   }
diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp
index 34fddf708a..97a127b318 100644
--- a/src/MOLECULE/improper_umbrella.cpp
+++ b/src/MOLECULE/improper_umbrella.cpp
@@ -25,10 +25,12 @@
 #include "domain.h"
 #include "force.h"
 #include "update.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 #define SMALL     0.001
@@ -269,7 +271,7 @@ void ImproperUmbrella::coeff(int narg, char **arg)
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     kw[i] = k_one;
-    w0[i] = w_one/180.0 * PI;
+    w0[i] = w_one/180.0 * MY_PI;
     if (w_one == 0) C[i] = 1.0;
     else C[i] = kw[i]/(pow(sin(w0[i]),2));
     setflag[i] = 1;
diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp
index eae9e4aa83..104cf83e6d 100644
--- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp
+++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp
@@ -26,11 +26,13 @@
 #include "neighbor.h"
 #include "neigh_request.h"
 #include "neigh_list.h"
+#include "domain.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
-#include "domain.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 #define CHUNK 8
@@ -44,8 +46,6 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp)
 
   no_virial_fdotr_compute = 1;
 
-  PI = 4.0*atan(1.0);
-
   nparams = maxparam = 0;
   params = NULL;
 
@@ -159,7 +159,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
 	  if (c < -1.0) c = -1.0;
 	  ac = acos(c);
 
-	  if (ac > pm->cut_angle && ac < (2.0*PI - pm->cut_angle)) {
+	  if (ac > pm->cut_angle && ac < (2.0*MY_PI - pm->cut_angle)) {
 	    s = sqrt(1.0 - c*c);
 	    if (s < SMALL) s = SMALL;
 
@@ -279,7 +279,7 @@ void PairHbondDreidingLJ::settings(int narg, char **arg)
   ap_global = force->inumeric(arg[0]);  
   cut_inner_global = force->numeric(arg[1]);
   cut_outer_global = force->numeric(arg[2]);
-  cut_angle_global = force->numeric(arg[3]) * PI/180.0;
+  cut_angle_global = force->numeric(arg[3]) * MY_PI/180.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -316,7 +316,7 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg)
   if (cut_inner_one>cut_outer_one)
     error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
   double cut_angle_one = cut_angle_global;
-  if (narg == 10) cut_angle_one = force->numeric(arg[9]) * PI/180.0;
+  if (narg == 10) cut_angle_one = force->numeric(arg[9]) * MY_PI/180.0;
   // grow params array if necessary
 
   if (nparams == maxparam) {
@@ -500,7 +500,7 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype,
     if (c < -1.0) c = -1.0;
     ac = acos(c);
 
-    if (ac < pm->cut_angle || ac > (2.0*PI - pm->cut_angle)) return 0.0;
+    if (ac < pm->cut_angle || ac > (2.0*MY_PI - pm->cut_angle)) return 0.0;
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;
 
diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.h b/src/MOLECULE/pair_hbond_dreiding_lj.h
index 387cd3f095..1304743b8d 100644
--- a/src/MOLECULE/pair_hbond_dreiding_lj.h
+++ b/src/MOLECULE/pair_hbond_dreiding_lj.h
@@ -38,7 +38,6 @@ class PairHbondDreidingLJ : public Pair {
  protected:
   double cut_inner_global,cut_outer_global,cut_angle_global;
   int ap_global;
-  double PI;
 
   struct Param {
     double epsilon,sigma;
diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp
index ebd3b27ef0..e7f9c1e7c5 100644
--- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp
+++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp
@@ -26,11 +26,13 @@
 #include "neighbor.h"
 #include "neigh_request.h"
 #include "neigh_list.h"
+#include "domain.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
-#include "domain.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 #define CHUNK 8
@@ -128,7 +130,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
 	  if (c < -1.0) c = -1.0;
 	  ac = acos(c);
 
-	  if (ac > pm->cut_angle && ac < (2.0*PI - pm->cut_angle)) {
+	  if (ac > pm->cut_angle && ac < (2.0*MY_PI - pm->cut_angle)) {
 	    s = sqrt(1.0 - c*c);
 	    if (s < SMALL) s = SMALL;
 
@@ -242,7 +244,7 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg)
   if (cut_inner_one>cut_outer_one)
     error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
   double cut_angle_one = cut_angle_global;
-  if (narg > 10) cut_angle_one = force->numeric(arg[10]) * PI/180.0;
+  if (narg > 10) cut_angle_one = force->numeric(arg[10]) * MY_PI/180.0;
 
   // grow params array if necessary
 
@@ -404,7 +406,7 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype,
     if (c < -1.0) c = -1.0;
     ac = acos(c);
 
-    if (ac < pm->cut_angle || ac > (2.0*PI - pm->cut_angle)) return 0.0;
+    if (ac < pm->cut_angle || ac > (2.0*MY_PI - pm->cut_angle)) return 0.0;
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;
 
diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp
index ec0d07d20d..411beb9c9d 100644
--- a/src/SRD/fix_srd.cpp
+++ b/src/SRD/fix_srd.cpp
@@ -34,10 +34,12 @@
 #include "fix_wall_srd.h"
 #include "random_mars.h"
 #include "random_park.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{SLIP,NOSLIP};
 enum{SPHERE,ELLIPSOID,WALL};
@@ -2087,8 +2089,6 @@ int FixSRD::update_srd(int i, double dt, double *xscoll, double *vsnew,
 
 void FixSRD::parameterize()
 {
-  double PI = 4.0*atan(1.0);
-
   // timesteps
 
   dt_big = update->dt;
@@ -2194,20 +2194,20 @@ void FixSRD::parameterize()
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & biggroupbit) {
 	if (radius && radius[i] > 0.0)
-	  volbig += 4.0/3.0*PI*radius[i]*radius[i]*radius[i];
+	  volbig += 4.0/3.0*MY_PI*radius[i]*radius[i]*radius[i];
 	else if (ellipsoid && ellipsoid[i] >= 0) {
 	  double *shape = ebonus[ellipsoid[i]].shape;
-	  volbig += 4.0/3.0*PI * shape[0]*shape[1]*shape[2];
+	  volbig += 4.0/3.0*MY_PI * shape[0]*shape[1]*shape[2];
 	}
       }
   } else {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & biggroupbit) {
 	if (radius && radius[i] > 0.0)
-	  volbig += PI*radius[i]*radius[i];
+	  volbig += MY_PI*radius[i]*radius[i];
 	else if (ellipsoid && ellipsoid[i] >= 0) {
 	  double *shape = ebonus[ellipsoid[i]].shape;
-	  volbig += PI*shape[0]*shape[1];
+	  volbig += MY_PI*shape[0]*shape[1];
 	}
       }
   }
diff --git a/src/USER-CG-CMM/angle_cg_cmm.cpp b/src/USER-CG-CMM/angle_cg_cmm.cpp
index 7d9800808a..070d65c5b8 100644
--- a/src/USER-CG-CMM/angle_cg_cmm.cpp
+++ b/src/USER-CG-CMM/angle_cg_cmm.cpp
@@ -24,10 +24,12 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 0.001
 
@@ -323,7 +325,7 @@ void AngleCGCMM::coeff(int narg, char **arg)
   for (int i = ilo; i <= ihi; i++) {
     k[i] = k_one;
     // convert theta0 from degrees to radians
-    theta0[i] = theta0_one/180.0 * PI;
+    theta0[i] = theta0_one/180.0 * MY_PI;
     epsilon[i] = epsilon_one;
     sigma[i] = sigma_one;
     rcut[i] = rcut_one; 
diff --git a/src/USER-CG-CMM/pair_cmm_common.cpp b/src/USER-CG-CMM/pair_cmm_common.cpp
index 81d142f568..4e1606d587 100644
--- a/src/USER-CG-CMM/pair_cmm_common.cpp
+++ b/src/USER-CG-CMM/pair_cmm_common.cpp
@@ -23,8 +23,10 @@
 #include "string.h"
 #include "ctype.h"
 #include "math.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 1.0e-6
 
@@ -302,15 +304,14 @@ double PairCMMCommon::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig2 = sigma[i][j]*sigma[i][j];
     double sig6 = sig2*sig2*sig2;
     double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
     double rc6 = rc3*rc3;
     double rc9 = rc3*rc6;
-    etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); 
-    ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); 
 #endif
   } 
diff --git a/src/USER-CUDA/fix_gravity_cuda.cpp b/src/USER-CUDA/fix_gravity_cuda.cpp
index 8c6d8488c8..36d7022fa3 100644
--- a/src/USER-CUDA/fix_gravity_cuda.cpp
+++ b/src/USER-CUDA/fix_gravity_cuda.cpp
@@ -30,12 +30,13 @@
 #include "update.h"
 #include "domain.h"
 #include "respa.h"
-#include "error.h"
 #include "cuda.h"
 #include "cuda_modify_flags.h"
-
+#include "math_const.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{CHUTE,SPHERICAL,GRADIENT,VECTOR};
 
@@ -79,8 +80,7 @@ FixGravityCuda::FixGravityCuda(LAMMPS *lmp, int narg, char **arg) :
     zdir = atof(arg[7]);
   } else error->all(FLERR,"Illegal fix gravity command");
 
-  double PI = 4.0*atan(1.0);
-  degree2rad = PI/180.0;
+  degree2rad = MY_PI/180.0;
 
   if (style == CHUTE || style == SPHERICAL || style == GRADIENT) {
     if (domain->dimension == 3) {
diff --git a/src/USER-CUDA/fix_shake_cuda.cpp b/src/USER-CUDA/fix_shake_cuda.cpp
index 219ac679bd..cbcb9acaef 100644
--- a/src/USER-CUDA/fix_shake_cuda.cpp
+++ b/src/USER-CUDA/fix_shake_cuda.cpp
@@ -34,8 +34,10 @@
 #include "error.h"
 #include "cuda.h"
 #include "cuda_modify_flags.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define BIG 1.0e20
 #define MASSDELTA 0.1
@@ -53,7 +55,6 @@ FixShakeCuda::FixShakeCuda(LAMMPS *lmp, int narg, char **arg) :
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
   neighbor_step=true;
-  PI = 4.0*atan(1.0);
 
   virial_flag = 1;
   create_attribute = 1;
@@ -2183,7 +2184,7 @@ void FixShakeCuda::stats()
       r3 = sqrt(delx*delx + dely*dely + delz*delz);
 
       angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
-      angle *= 180.0/PI;
+      angle *= 180.0/MY_PI;
       m = shake_type[i][2];
       a_count[m]++;
       a_ave[m] += angle;
diff --git a/src/USER-CUDA/fix_shake_cuda.h b/src/USER-CUDA/fix_shake_cuda.h
index 18ea64f983..b0ebe584d7 100644
--- a/src/USER-CUDA/fix_shake_cuda.h
+++ b/src/USER-CUDA/fix_shake_cuda.h
@@ -53,7 +53,6 @@ class FixShakeCuda : public Fix {
  private:
   class Cuda *cuda;
   int me,nprocs;
-  double PI;
   double tolerance;                      // SHAKE tolerance
   int max_iter;                          // max # of SHAKE iterations
   int output_every;                      // SHAKE stat output every so often
diff --git a/src/USER-CUDA/pppm_cuda.cpp b/src/USER-CUDA/pppm_cuda.cpp
index 8167eb1d72..5fd8406fbe 100644
--- a/src/USER-CUDA/pppm_cuda.cpp
+++ b/src/USER-CUDA/pppm_cuda.cpp
@@ -61,8 +61,10 @@
 #include "cuda_wrapper_cu.h"
 #include "pppm_cuda_cu.h"
 #include "cuda.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define MAXORDER 7
 #define OFFSET 4096
@@ -109,7 +111,6 @@ PPPMCuda::PPPMCuda(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, (narg==2?1:nar
   if(narg>1)
   precisionmodify=arg[1][0];
   else precisionmodify='=';
-  PI = 4.0*atan(1.0);
   
   nfactors = 3;
   factors = new int[nfactors];
@@ -648,9 +649,9 @@ void PPPMCuda::setup()
 
   delvolinv = delxinv*delyinv*delzinv;
 
-  double unitkx = (2.0*PI/xprd);
-  double unitky = (2.0*PI/yprd);
-  double unitkz = (2.0*PI/zprd_slab);
+  double unitkx = (2.0*MY_PI/xprd);
+  double unitky = (2.0*MY_PI/yprd);
+  double unitkz = (2.0*MY_PI/zprd_slab);
 
   // fkx,fky,fkz for my FFT grid pts
   Cuda_PPPM_Setup_fkxyz_vg(nx_pppm, ny_pppm,nz_pppm,unitkx,unitky,unitkz,g_ewald);
@@ -747,11 +748,11 @@ double sqk;
   double sum1,dot1,dot2;
   double numerator,denominator;
 
-  int nbx = static_cast<int> ((g_ewald*xprd/(PI*nx_pppm)) * 
+  int nbx = static_cast<int> ((g_ewald*xprd/(MY_PI*nx_pppm)) * 
 			      pow(-log(EPS_HOC),0.25));
-  int nby = static_cast<int> ((g_ewald*yprd/(PI*ny_pppm)) * 
+  int nby = static_cast<int> ((g_ewald*yprd/(MY_PI*ny_pppm)) * 
 			      pow(-log(EPS_HOC),0.25));
-  int nbz = static_cast<int> ((g_ewald*zprd_slab/(PI*nz_pppm)) * 
+  int nbz = static_cast<int> ((g_ewald*zprd_slab/(MY_PI*nz_pppm)) * 
 			      pow(-log(EPS_HOC),0.25));
   Cuda_PPPM_setup_greensfn(nx_pppm,ny_pppm,nz_pppm,unitkx,unitky,unitkz,g_ewald,
 nbx,nby,nbz,xprd,yprd,zprd_slab);
@@ -967,7 +968,7 @@ else
    
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/1.772453851 +
-      0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
+      MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e;
   }
 
@@ -1728,11 +1729,11 @@ void PPPMCuda::slabcorr(int eflag)
 
   // compute corrections
   
-  double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume;
+  double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
   
   if (eflag) energy += qqrd2e*scale * e_slabcorr;
   
-  double ffact = -4.0*PI*dipole_all/volume; 
+  double ffact = -4.0*MY_PI*dipole_all/volume; 
  
   cuda_slabcorr_force(&cuda->shared_data,ffact);
 }
diff --git a/src/USER-EFF/atom_vec_electron.h b/src/USER-EFF/atom_vec_electron.h
index b9372f62a6..6c527cacdd 100644
--- a/src/USER-EFF/atom_vec_electron.h
+++ b/src/USER-EFF/atom_vec_electron.h
@@ -60,7 +60,6 @@ class AtomVecElectron : public AtomVec {
   bigint memory_usage();
   
  private:
-  double PI;
   int *tag,*type,*mask,*image;
   double **x,**v,**f;
   int *spin;
diff --git a/src/USER-OMP/dihedral_helix_omp.cpp b/src/USER-OMP/dihedral_helix_omp.cpp
index a3ca969ef3..4ec701a0cb 100644
--- a/src/USER-OMP/dihedral_helix_omp.cpp
+++ b/src/USER-OMP/dihedral_helix_omp.cpp
@@ -25,9 +25,11 @@
 #include "domain.h"
 #include "force.h"
 #include "update.h"
+#include "math_const.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 #define SMALL     0.001
@@ -211,10 +213,10 @@ void DihedralHelixOMP::eval(double **f, int nfrom, int nto, int tid)
     siinv = 1.0/si;
 
     pd = -aphi[type] + 3.0*bphi[type]*sin(3.0*phi)*siinv +
-      cphi[type]*sin(phi + 0.25*PI)*siinv;
+      cphi[type]*sin(phi + MY_PI4)*siinv;
 
     if (EFLAG) edihedral = aphi[type]*(1.0 - c) + bphi[type]*(1.0 + cos(3.0*phi)) +
-      cphi[type]*(1.0 + cos(phi + 0.25*PI));
+		 cphi[type]*(1.0 + cos(phi + MY_PI4));
 ;
 
     a = pd;
diff --git a/src/USER-OMP/pair_soft_omp.cpp b/src/USER-OMP/pair_soft_omp.cpp
index 7667efa983..9f9673a28b 100644
--- a/src/USER-OMP/pair_soft_omp.cpp
+++ b/src/USER-OMP/pair_soft_omp.cpp
@@ -19,8 +19,10 @@
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define SMALL 1.0e-4
 
@@ -122,7 +124,7 @@ void PairSoftOMP::eval(double **f, int iifrom, int iito, int tid)
 
       if (rsq < cutsq[itype][jtype]) {
 	r = sqrt(rsq);
-	arg = PI/cut[itype][jtype];
+	arg = MY_PI/cut[itype][jtype];
 	if (r > SMALL) fpair = factor_lj * prefactor[itype][jtype] * 
 		       sin(arg*r) * arg/r;
 	else fpair = 0.0;
diff --git a/src/angle.cpp b/src/angle.cpp
index fbed9116a9..58ba679948 100644
--- a/src/angle.cpp
+++ b/src/angle.cpp
@@ -15,10 +15,12 @@
 #include "angle.h"
 #include "atom.h"
 #include "force.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -27,8 +29,6 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
   energy = 0.0;
 
   allocated = 0;
-  PI = 4.0*atan(1.0);
-  THIRD = 1.0/3.0;
 
   maxeatom = maxvatom = 0;
   eatom = NULL;
diff --git a/src/angle.h b/src/angle.h
index 99d8dc105f..7f25f38289 100644
--- a/src/angle.h
+++ b/src/angle.h
@@ -40,8 +40,6 @@ class Angle : protected Pointers {
   virtual double memory_usage();
 
  protected:
-  double PI,THIRD;
-
   int evflag;
   int eflag_either,eflag_global,eflag_atom;
   int vflag_either,vflag_global,vflag_atom;
diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp
index 47c9ef97f6..4baff5b602 100755
--- a/src/atom_vec_ellipsoid.cpp
+++ b/src/atom_vec_ellipsoid.cpp
@@ -25,10 +25,12 @@
 #include "domain.h"
 #include "modify.h"
 #include "fix.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define DELTA 10000
 #define DELTA_BONUS 10000
@@ -53,8 +55,6 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) :
   atom->ellipsoid_flag = 1;
   atom->rmass_flag = atom->angmom_flag = atom->torque_flag = 1;
 
-  PI = 4.0*atan(1.0);
-
   nlocal_bonus = nghost_bonus = nmax_bonus = 0;
   bonus = NULL;
 }
@@ -1200,7 +1200,7 @@ void AtomVecEllipsoid::data_atom_bonus(int m, char **values)
   // reset ellipsoid mass
   // previously stored density in rmass
 
-  rmass[m] *= 4.0*PI/3.0 * shape[0]*shape[1]*shape[2];
+  rmass[m] *= 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2];
 
   bonus[nlocal_bonus].ilocal = m;
   ellipsoid[m] = nlocal_bonus++;
diff --git a/src/atom_vec_ellipsoid.h b/src/atom_vec_ellipsoid.h
index 537f1eaab7..0020eb5f28 100755
--- a/src/atom_vec_ellipsoid.h
+++ b/src/atom_vec_ellipsoid.h
@@ -76,7 +76,6 @@ class AtomVecEllipsoid : public AtomVec {
   void set_shape(int, double, double, double);
 
  private:
-  double PI;
   int *tag,*type,*mask,*image;
   double **x,**v,**f;
   double *rmass;
diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp
index a3e2c34772..96b2d766d8 100644
--- a/src/atom_vec_sphere.cpp
+++ b/src/atom_vec_sphere.cpp
@@ -23,10 +23,12 @@
 #include "force.h"
 #include "fix.h"
 #include "fix_adapt.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define DELTA 10000
 
@@ -50,8 +52,6 @@ AtomVecSphere::AtomVecSphere(LAMMPS *lmp, int narg, char **arg) :
   atom->sphere_flag = 1;
   atom->radius_flag = atom->rmass_flag = atom->omega_flag = 
     atom->torque_flag = 1;
-
-  PI = 4.0*atan(1.0);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -929,7 +929,7 @@ void AtomVecSphere::create_atom(int itype, double *coord)
   v[nlocal][2] = 0.0;
 
   radius[nlocal] = 0.5;
-  rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal];
+  rmass[nlocal] = 4.0*MY_PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal];
   omega[nlocal][0] = 0.0;
   omega[nlocal][1] = 0.0;
   omega[nlocal][2] = 0.0;
@@ -965,7 +965,7 @@ void AtomVecSphere::data_atom(double *coord, int imagetmp, char **values)
 
   if (radius[nlocal] == 0.0) rmass[nlocal] = density;
   else 
-    rmass[nlocal] = 4.0*PI/3.0 *
+    rmass[nlocal] = 4.0*MY_PI/3.0 *
       radius[nlocal]*radius[nlocal]*radius[nlocal] * density;
 
   x[nlocal][0] = coord[0];
@@ -1002,7 +1002,7 @@ int AtomVecSphere::data_atom_hybrid(int nlocal, char **values)
 
   if (radius[nlocal] == 0.0) rmass[nlocal] = density;
   else 
-    rmass[nlocal] = 4.0*PI/3.0 *
+    rmass[nlocal] = 4.0*MY_PI/3.0 *
       radius[nlocal]*radius[nlocal]*radius[nlocal] * density;
 
   return 2;
diff --git a/src/atom_vec_sphere.h b/src/atom_vec_sphere.h
index 8157029615..78357b7d1a 100644
--- a/src/atom_vec_sphere.h
+++ b/src/atom_vec_sphere.h
@@ -61,7 +61,6 @@ class AtomVecSphere : public AtomVec {
   bigint memory_usage();
 
  private:
-  double PI;
   int *tag,*type,*mask,*image;
   double **x,**v,**f;
   double *radius,*density,*rmass;
diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp
index a126348d31..d19b7e71aa 100644
--- a/src/compute_angle_local.cpp
+++ b/src/compute_angle_local.cpp
@@ -20,10 +20,12 @@
 #include "domain.h"
 #include "force.h"
 #include "angle.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define DELTA 10000
 
@@ -133,7 +135,6 @@ int ComputeAngleLocal::compute_angles(int flag)
   }
 
   Angle *angle = force->angle;
-  double PI = 4.0*atan(1.0);
 
   m = n = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
@@ -170,7 +171,7 @@ int ComputeAngleLocal::compute_angles(int flag)
 	  c /= r1*r2;
 	  if (c > 1.0) c = 1.0;
 	  if (c < -1.0) c = -1.0;
-	  tbuf[n] = 180.0*acos(c)/PI;
+	  tbuf[n] = 180.0*acos(c)/MY_PI;
 	}
 
 	if (eflag >= 0) {
diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp
index 6215346717..10cac6f607 100644
--- a/src/compute_dihedral_local.cpp
+++ b/src/compute_dihedral_local.cpp
@@ -20,10 +20,12 @@
 #include "domain.h"
 #include "force.h"
 #include "dihedral.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define DELTA 10000
 #define SMALL 0.001
@@ -128,8 +130,6 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
     }
   }
 
-  double PI = 4.0*atan(1.0);
-
   m = n = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
     if (!(mask[atom2] & groupbit)) continue;
@@ -190,7 +190,7 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
 
 	  if (c > 1.0) c = 1.0;
 	  if (c < -1.0) c = -1.0;
-	  pbuf[n] = 180.0*atan2(s,c)/PI;
+	  pbuf[n] = 180.0*atan2(s,c)/MY_PI;
 	}
 	n += nvalues;
       }
diff --git a/src/compute_improper_local.cpp b/src/compute_improper_local.cpp
index 7162277c39..3285f2561c 100644
--- a/src/compute_improper_local.cpp
+++ b/src/compute_improper_local.cpp
@@ -20,10 +20,12 @@
 #include "domain.h"
 #include "force.h"
 #include "improper.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define DELTA 10000
 
@@ -129,8 +131,6 @@ int ComputeImproperLocal::compute_impropers(int flag)
     }
   }
 
-  double PI = 4.0*atan(1.0);
-
   m = n = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
     if (!(mask[atom2] & groupbit)) continue;
@@ -188,7 +188,7 @@ int ComputeImproperLocal::compute_impropers(int flag)
 
 	  if (c > 1.0) c = 1.0;
 	  if (c < -1.0) c = -1.0;
-	  cbuf[n] = 180.0*acos(c)/PI;
+	  cbuf[n] = 180.0*acos(c)/MY_PI;
 	}
 	n += nvalues;
       }
diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp
index 9dbe3aded0..31776f9473 100644
--- a/src/compute_rdf.cpp
+++ b/src/compute_rdf.cpp
@@ -28,10 +28,12 @@
 #include "neigh_request.h"
 #include "neigh_list.h"
 #include "group.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -267,10 +269,9 @@ void ComputeRDF::compute_array()
   //   assuming J atoms are at uniform density
 
   double constant,nideal,gr,ncoord,rlower,rupper;
-  double PI = 4.0*atan(1.0);
 
   if (domain->dimension == 3) {
-    constant = 4.0*PI / (3.0*domain->xprd*domain->yprd*domain->zprd);
+    constant = 4.0*MY_PI / (3.0*domain->xprd*domain->yprd*domain->zprd);
 
     for (m = 0; m < npairs; m++) {
       ncoord = 0.0;
@@ -289,7 +290,7 @@ void ComputeRDF::compute_array()
     }
 
   } else {
-    constant = PI / (domain->xprd*domain->yprd);
+    constant = MY_PI / (domain->xprd*domain->yprd);
 
     for (m = 0; m < npairs; m++) {
       ncoord = 0.0;
diff --git a/src/dihedral.cpp b/src/dihedral.cpp
index a3edadc0c0..974df81b14 100644
--- a/src/dihedral.cpp
+++ b/src/dihedral.cpp
@@ -32,7 +32,6 @@ Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp)
   energy = 0.0;
 
   allocated = 0;
-  PI = 4.0*atan(1.0);
 
   maxeatom = maxvatom = 0;
   eatom = NULL;
diff --git a/src/dihedral.h b/src/dihedral.h
index c41cedad03..b41272b7ba 100644
--- a/src/dihedral.h
+++ b/src/dihedral.h
@@ -41,8 +41,6 @@ class Dihedral : protected Pointers {
   virtual double memory_usage();
 
  protected:
-  double PI;
-
   int evflag;
   int eflag_either,eflag_global,eflag_atom;
   int vflag_either,vflag_global,vflag_atom;
diff --git a/src/dump_image.cpp b/src/dump_image.cpp
index 25dc36109f..16f677327e 100644
--- a/src/dump_image.cpp
+++ b/src/dump_image.cpp
@@ -30,6 +30,7 @@
 #include "input.h"
 #include "variable.h"
 #include "random_mars.h"
+#include "math_const.h"
 #include "error.h"
 #include "memory.h"
 
@@ -38,6 +39,7 @@
 #endif
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define NCOLORS 140
 #define NELEMENTS 109
@@ -57,8 +59,6 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 {
   if (binary || multiproc) error->all(FLERR,"Invalid dump image filename");
 
-  PI = 4.0*atan(1.0);
-
   // set filetype based on filename suffix
 
   int n = strlen(filename);
@@ -95,8 +95,8 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
     bdiamvalue = 0.5;
   }
   width = height = 512;
-  theta = 60.0 * PI/180.0;
-  phi = 30.0 * PI/180.0;
+  theta = 60.0 * MY_PI/180.0;
+  phi = 30.0 * MY_PI/180.0;
   thetastr = phistr = NULL;
   cflag = STATIC;
   cx = cy = cz = 0.5;
@@ -171,7 +171,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 	theta = atof(arg[iarg+1]);
 	if (theta < 0.0 || theta > 180.0)
 	  error->all(FLERR,"Invalid dump image theta value");
-	theta *= PI/180.0;
+	theta *= MY_PI/180.0;
       }
       if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
 	int n = strlen(&arg[iarg+2][2]) + 1;
@@ -179,7 +179,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 	strcpy(phistr,&arg[iarg+2][2]);
       } else {
 	phi = atof(arg[iarg+2]);
-	phi *= PI/180.0;
+	phi *= MY_PI/180.0;
       }
       iarg += 3;
 
@@ -358,25 +358,25 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 
   // static parameters
 
-  FOV = PI/6.0;              // 30 degrees
+  FOV = MY_PI/6.0;              // 30 degrees
   ambientColor[0] = 0.0;
   ambientColor[1] = 0.0;
   ambientColor[2] = 0.0;
 
-  keyLightPhi = -PI/4.0;     // -45 degrees
-  keyLightTheta = PI/6.0;    // 30 degrees
+  keyLightPhi = -MY_PI4;     // -45 degrees
+  keyLightTheta = MY_PI/6.0;    // 30 degrees
   keyLightColor[0] = 0.9;
   keyLightColor[1] = 0.9;
   keyLightColor[2] = 0.9;
 
-  fillLightPhi = PI/6.0;     // 30 degrees
+  fillLightPhi = MY_PI/6.0;     // 30 degrees
   fillLightTheta = 0; 
   fillLightColor[0] = 0.45;
   fillLightColor[1] = 0.45;
   fillLightColor[2] = 0.45;
 
-  backLightPhi = PI;         // 180 degrees
-  backLightTheta = PI/12.0;  // 15 degrees
+  backLightPhi = MY_PI;         // 180 degrees
+  backLightTheta = MY_PI/12.0;  // 15 degrees
   backLightColor[0] = 0.9;
   backLightColor[1] = 0.9;
   backLightColor[2] = 0.9;
@@ -676,11 +676,11 @@ void DumpImage::view_params()
     theta = input->variable->compute_equal(thetavar);
     if (theta < 0.0 || theta > 180.0)
       error->all(FLERR,"Invalid dump image theta value");
-    theta *= PI/180.0;
+    theta *= MY_PI/180.0;
   }
   if (phistr) {
     phi = input->variable->compute_equal(phivar);
-    phi *= PI/180.0;
+    phi *= MY_PI/180.0;
   }
 
   camDir[0] = sin(theta)*cos(phi);
@@ -764,7 +764,7 @@ void DumpImage::view_params()
   if (ssao) {
     SSAORadius = maxdel * 0.05 * ssaoint;
     SSAOSamples = static_cast<int> (8.0 + 32.0*ssaoint);
-    SSAOJitter = PI / 12;
+    SSAOJitter = MY_PI / 12;
     ambientColor[0] = 0.5;
     ambientColor[1] = 0.5;
     ambientColor[2] = 0.5;
@@ -1358,7 +1358,7 @@ void DumpImage::compute_SSAO()
 {
   // used for rasterizing the spheres
 
-  double delTheta = 2.0*PI / SSAOSamples;
+  double delTheta = 2.0*MY_PI / SSAOSamples;
 
   // typical neighborhood value for shading
 
diff --git a/src/dump_image.h b/src/dump_image.h
index 866fc32e89..844ba1c112 100644
--- a/src/dump_image.h
+++ b/src/dump_image.h
@@ -69,8 +69,6 @@ class DumpImage : public DumpCustom {
 
   // constant view params
 
-  double PI;
-
   double FOV;
   double ambientColor[3];
 
diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
index fd1a256b8a..59af597bc3 100644
--- a/src/fix_adapt.cpp
+++ b/src/fix_adapt.cpp
@@ -24,10 +24,12 @@
 #include "kspace.h"
 #include "input.h"
 #include "variable.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{PAIR,KSPACE,ATOM};
 enum{DIAMETER};
@@ -309,7 +311,6 @@ void FixAdapt::change_settings()
       if (ad->aparam == DIAMETER) {
 	int mflag = 0;
 	if (atom->rmass_flag) mflag = 1;
-	double PI = 4.0*atan(1.0);
 	double density;
 
 	double *radius = atom->radius;
@@ -324,9 +325,11 @@ void FixAdapt::change_settings()
 	} else {
 	  for (i = 0; i < nlocal; i++)
 	    if (mask[i] & groupbit) {
-	      density = rmass[i] / (4.0*PI/3.0 * radius[i]*radius[i]*radius[i]);
+	      density = rmass[i] / (4.0*MY_PI/3.0 * 
+				    radius[i]*radius[i]*radius[i]);
 	      radius[i] = 0.5*value;
-	      rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density;
+	      rmass[i] = 4.0*MY_PI/3.0 * 
+		radius[i]*radius[i]*radius[i] * density;
 	    }
 	}
       }
diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp
index 039bba40b2..097683ad23 100644
--- a/src/fix_deform.cpp
+++ b/src/fix_deform.cpp
@@ -27,10 +27,12 @@
 #include "lattice.h"
 #include "force.h"
 #include "modify.h"
+#include "math_const.h"
 #include "kspace.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE};
 enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE};
@@ -305,7 +307,7 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   if (force_reneighbor) irregular = new Irregular(lmp);
   else irregular = NULL;
 
-  TWOPI = 8.0*atan(1.0);
+  TWOPI = 2.0*MY_PI;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp
index 09a2e2c7f1..de71e6bd79 100644
--- a/src/fix_gravity.cpp
+++ b/src/fix_gravity.cpp
@@ -20,9 +20,11 @@
 #include "update.h"
 #include "domain.h"
 #include "respa.h"
+#include "math_const.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{CHUTE,SPHERICAL,GRADIENT,VECTOR};
 
@@ -65,8 +67,7 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
     zdir = atof(arg[7]);
   } else error->all(FLERR,"Illegal fix gravity command");
 
-  double PI = 4.0*atan(1.0);
-  degree2rad = PI/180.0;
+  degree2rad = MY_PI/180.0;
 
   if (style == CHUTE || style == SPHERICAL || style == GRADIENT) {
     if (domain->dimension == 3) {
diff --git a/src/fix_move.cpp b/src/fix_move.cpp
index afd087772e..74657fa86e 100644
--- a/src/fix_move.cpp
+++ b/src/fix_move.cpp
@@ -26,10 +26,12 @@
 #include "respa.h"
 #include "input.h"
 #include "variable.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
 enum{EQUAL,ATOM};
@@ -216,10 +218,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
 
   // set omega_rotate from period
 
-  if (mstyle == WIGGLE || mstyle == ROTATE) {
-    double PI = 4.0 * atan(1.0);
-    omega_rotate = 2.0*PI / period;
-  }
+  if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = 2.0*MY_PI / period;
 
   // runit = unit vector along rotation axis
 
diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp
index d0269acddb..109373add7 100644
--- a/src/fix_orient_fcc.cpp
+++ b/src/fix_orient_fcc.cpp
@@ -28,10 +28,12 @@
 #include "neigh_request.h"
 #include "comm.h"
 #include "output.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define BIG 1000000000
 
@@ -77,7 +79,6 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
 
   // initializations
 
-  PI = 4.0*atan(1.0);
   half_fcc_nn = 6;
   use_xismooth = false;
   double xicutoff = 1.57;
@@ -344,10 +345,10 @@ void FixOrientFCC::post_force(int vflag)
       edelta = Vxi;
       order[i][1] = 1.0;
     } else {
-      omega = (0.5*PI)*(xi_total-xi0) / (xi1-xi0);
-      nbr[i].duxi = PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0));
+      omega = MY_PI2*(xi_total-xi0) / (xi1-xi0);
+      nbr[i].duxi = MY_PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0));
       edelta = Vxi*(1 - cos(2.0*omega)) / 2.0;
-      order[i][1] = omega / (0.5*PI);
+      order[i][1] = omega / MY_PI2;
     }
     added_energy += edelta;
   }
diff --git a/src/fix_orient_fcc.h b/src/fix_orient_fcc.h
index f401d30855..b83753619b 100644
--- a/src/fix_orient_fcc.h
+++ b/src/fix_orient_fcc.h
@@ -58,7 +58,6 @@ class FixOrientFCC : public Fix {
 
  private:
   int me;
-  double PI;
   int nlevels_respa;
 
   int direction_of_motion;         // 1 = center shrinks, 0 = center grows
diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp
index 43f7f61087..f3048cc18d 100644
--- a/src/fix_restrain.cpp
+++ b/src/fix_restrain.cpp
@@ -26,10 +26,12 @@
 #include "comm.h"
 #include "respa.h"
 #include "input.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{DIHEDRAL};
 
@@ -83,13 +85,12 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
   // special treatment for dihedral restraints
 
   if (rstyle == DIHEDRAL) {
-    double PI = 4.0*atan(1.0);
     cos_shift = (double *)
       memory->smalloc(n_bonds * sizeof(double), "restrain:cos_shift");
     sin_shift = (double *)
       memory->smalloc(n_bonds * sizeof(double), "restrain:sin_shift");
     for (ibond = 0; ibond < n_bonds; ibond++) {
-      double my_arg = PI * (180.0 + target[ibond]) / 180.0;
+      double my_arg = MY_PI * (180.0 + target[ibond]) / 180.0;
       cos_shift[ibond] = cos(my_arg);
       sin_shift[ibond] = sin(my_arg);
     }
diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp
index e00e84e424..5cdce3ad25 100644
--- a/src/fix_shake.cpp
+++ b/src/fix_shake.cpp
@@ -30,10 +30,12 @@
 #include "comm.h"
 #include "group.h"
 #include "fix_respa.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define BIG 1.0e20
 #define MASSDELTA 0.1
@@ -46,8 +48,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
 
-  PI = 4.0*atan(1.0);
-
   virial_flag = 1;
   create_attribute = 1;
 
@@ -2114,7 +2114,7 @@ void FixShake::stats()
       r3 = sqrt(delx*delx + dely*dely + delz*delz);
 
       angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
-      angle *= 180.0/PI;
+      angle *= 180.0/MY_PI;
       m = shake_type[i][2];
       a_count[m]++;
       a_ave[m] += angle;
diff --git a/src/fix_shake.h b/src/fix_shake.h
index eaf39837cb..b0aa5c63d2 100644
--- a/src/fix_shake.h
+++ b/src/fix_shake.h
@@ -49,7 +49,6 @@ class FixShake : public Fix {
 
  private:
   int me,nprocs;
-  double PI;
   double tolerance;                      // SHAKE tolerance
   int max_iter;                          // max # of SHAKE iterations
   int output_every;                      // SHAKE stat output every so often
diff --git a/src/improper.cpp b/src/improper.cpp
index 3c65e8c89c..8f5ab4d6e3 100644
--- a/src/improper.cpp
+++ b/src/improper.cpp
@@ -27,7 +27,6 @@ Improper::Improper(LAMMPS *lmp) : Pointers(lmp)
   energy = 0.0;
 
   allocated = 0;
-  PI = 4.0*atan(1.0);
 
   maxeatom = maxvatom = 0;
   eatom = NULL;
diff --git a/src/improper.h b/src/improper.h
index 711f124c20..d8310b2475 100644
--- a/src/improper.h
+++ b/src/improper.h
@@ -38,8 +38,6 @@ class Improper : protected Pointers {
   virtual double memory_usage();
 
  protected:
-  double PI;
-
   int evflag;
   int eflag_either,eflag_global,eflag_atom;
   int vflag_either,vflag_global,vflag_atom;
diff --git a/src/pair_born.cpp b/src/pair_born.cpp
index 406b95a0cc..271b7b1f7a 100644
--- a/src/pair_born.cpp
+++ b/src/pair_born.cpp
@@ -24,10 +24,12 @@
 #include "comm.h"
 #include "force.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -269,7 +271,6 @@ double PairBorn::init_one(int i, int j)
      } 
      MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
 
-     double PI = 4.0*atan(1.0);
      double rho1 = rho[i][j];
      double rho2 = rho1*rho1;
      double rho3 = rho2*rho1;
@@ -277,11 +278,11 @@ double PairBorn::init_one(int i, int j)
      double rc2 = rc*rc;
      double rc3 = rc2*rc;
      double rc5 = rc3*rc2;
-     etail_ij = 2.0*PI*all[0]*all[1] * 
+     etail_ij = 2.0*MY_PI*all[0]*all[1] * 
        (a[i][j]*exp((sigma[i][j]-rc)/rho1)*rho1* 
 	(rc2 + 2.0*rho1*rc + 2.0*rho2) - 
 	c[i][j]/(3.0*rc3) + d[i][j]/(5.0*rc5));
-     ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1] * 
+     ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1] * 
        (-a[i][j]*exp((sigma[i][j]-rc)/rho1) * 
 	(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 
 	2.0*c[i][j]/rc3 - 8.0*d[i][j]/(5.0*rc5)); 
diff --git a/src/pair_buck.cpp b/src/pair_buck.cpp
index 0490f80eea..1a3464e27e 100644
--- a/src/pair_buck.cpp
+++ b/src/pair_buck.cpp
@@ -20,10 +20,12 @@
 #include "comm.h"
 #include "force.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -249,17 +251,16 @@ double PairBuck::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
 
-    double PI = 4.0*atan(1.0);
     double rho1 = rho[i][j];
     double rho2 = rho1*rho1;
     double rho3 = rho2*rho1;
     double rc = cut[i][j];
     double rc2 = rc*rc;
     double rc3 = rc2*rc;
-    etail_ij = 2.0*PI*all[0]*all[1]*
+    etail_ij = 2.0*MY_PI*all[0]*all[1]*
       (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) - 
        c[i][j]/(3.0*rc3));
-    ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1]*
+    ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1]*
       (-a[i][j]*exp(-rc/rho1)*
        (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3);
   }
diff --git a/src/pair_buck_coul_cut.cpp b/src/pair_buck_coul_cut.cpp
index 51c179b889..a2f0784a3d 100644
--- a/src/pair_buck_coul_cut.cpp
+++ b/src/pair_buck_coul_cut.cpp
@@ -24,10 +24,12 @@
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -304,17 +306,16 @@ double PairBuckCoulCut::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
 
-    double PI = 4.0*atan(1.0);
     double rho1 = rho[i][j];
     double rho2 = rho1*rho1;
     double rho3 = rho2*rho1;
     double rc = cut_lj[i][j];
     double rc2 = rc*rc;
     double rc3 = rc2*rc;
-    etail_ij = 2.0*PI*all[0]*all[1]*
+    etail_ij = 2.0*MY_PI*all[0]*all[1]*
       (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) - 
        c[i][j]/(3.0*rc3));
-    ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1]*
+    ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1]*
       (-a[i][j]*exp(-rc/rho1)*
        (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3);
   }
diff --git a/src/pair_lj96_cut.cpp b/src/pair_lj96_cut.cpp
index 76f808dd03..8eb9a3d9d5 100644
--- a/src/pair_lj96_cut.cpp
+++ b/src/pair_lj96_cut.cpp
@@ -29,10 +29,12 @@
 #include "update.h"
 #include "integrate.h"
 #include "respa.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -589,15 +591,14 @@ double PairLJ96Cut::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
     double sig6 = sig3*sig3;
     double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
     double rc6 = rc3*rc3;
     
-    etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig3 - 2.0*rc3) / (6.0*rc6); 
-    ptail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (3.0*sig3 - 4.0*rc3) / (6.0*rc6); 
   } 
 
diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp
index 671d6420cf..44d86990df 100644
--- a/src/pair_lj_cut.cpp
+++ b/src/pair_lj_cut.cpp
@@ -29,10 +29,12 @@
 #include "update.h"
 #include "integrate.h"
 #include "respa.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -583,15 +585,14 @@ double PairLJCut::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig2 = sigma[i][j]*sigma[i][j];
     double sig6 = sig2*sig2*sig2;
     double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
     double rc6 = rc3*rc3;
     double rc9 = rc3*rc6;
-    etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); 
-    ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); 
   } 
 
diff --git a/src/pair_lj_cut_coul_cut.cpp b/src/pair_lj_cut_coul_cut.cpp
index 1d311d4b72..c6973b2580 100644
--- a/src/pair_lj_cut_coul_cut.cpp
+++ b/src/pair_lj_cut_coul_cut.cpp
@@ -21,10 +21,12 @@
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -300,15 +302,14 @@ double PairLJCutCoulCut::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
         
-    double PI = 4.0*atan(1.0);
     double sig2 = sigma[i][j]*sigma[i][j];
     double sig6 = sig2*sig2*sig2;
     double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
     double rc6 = rc3*rc3;
     double rc9 = rc3*rc6;
-    etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); 
-    ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); 
   } 
 
diff --git a/src/pair_lj_expand.cpp b/src/pair_lj_expand.cpp
index 57110ac81c..10e69ba977 100644
--- a/src/pair_lj_expand.cpp
+++ b/src/pair_lj_expand.cpp
@@ -19,10 +19,12 @@
 #include "comm.h"
 #include "force.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -259,7 +261,6 @@ double PairLJExpand::init_one(int i, int j)
     }
     MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
 
-    double PI = 4.0*atan(1.0);
     double sig2 = sigma[i][j]*sigma[i][j];
     double sig6 = sig2*sig2*sig2;
     double shiftcut = shift[i][j] - cut[i][j];
@@ -273,11 +274,11 @@ double PairLJExpand::init_one(int i, int j)
     double rc12 = rc11*shiftcut;
     double shift2 = shift[i][j]*shift[i][j];
     double shift3 = shift2*shift[i][j];
-    etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6*((-1.0/(9.0*rc9) + shift[i][j]/(5.0*rc10) -
              shift2/(11.0*rc11))*sig6 + 
 	    1.0/(3.0*rc3) - shift[i][j]/(2.0*rc4) + shift2/(5.0*rc5));
-    ptail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * 
+    ptail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * 
       sig6* ((-4.0/(3.0*rc9) + 18.0*shift[i][j]/(5.0*rc10) -
 	     36.0*shift2/(11.0*rc11) + shift3/rc12)*sig6 + 
 	     2.0/rc3 - 9.0*shift[i][j]/(2.0*rc4) + 
diff --git a/src/pair_soft.cpp b/src/pair_soft.cpp
index ce25f4a782..63f0b72b05 100644
--- a/src/pair_soft.cpp
+++ b/src/pair_soft.cpp
@@ -21,17 +21,16 @@
 #include "force.h"
 #include "update.h"
 #include "neigh_list.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-PairSoft::PairSoft(LAMMPS *lmp) : Pair(lmp)
-{
-  PI = 4.0*atan(1.0);
-}
+PairSoft::PairSoft(LAMMPS *lmp) : Pair(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
@@ -95,9 +94,9 @@ void PairSoft::compute(int eflag, int vflag)
 
       if (rsq < cutsq[itype][jtype]) {
 	r = sqrt(rsq);
-	arg = PI*r/cut[itype][jtype];
+	arg = MY_PI*r/cut[itype][jtype];
 	if (r > 0.0) fpair = factor_lj * prefactor[itype][jtype] * 
-		       sin(arg) * PI/cut[itype][jtype]/r;
+		       sin(arg) * MY_PI/cut[itype][jtype]/r;
 	else fpair = 0.0;
 
 	f[i][0] += delx*fpair;
@@ -290,9 +289,9 @@ double PairSoft::single(int i, int j, int itype, int jtype, double rsq,
   double r,arg,philj;
 
   r = sqrt(rsq);
-  arg = PI*r/cut[itype][jtype];
+  arg = MY_PI*r/cut[itype][jtype];
   fforce = factor_lj * prefactor[itype][jtype] * 
-    sin(arg) * PI/cut[itype][jtype]/r;
+    sin(arg) * MY_PI/cut[itype][jtype]/r;
   
   philj = prefactor[itype][jtype] * (1.0+cos(arg));
   return factor_lj*philj;
diff --git a/src/pair_soft.h b/src/pair_soft.h
index 1558f899bc..5cb8950659 100644
--- a/src/pair_soft.h
+++ b/src/pair_soft.h
@@ -43,7 +43,6 @@ class PairSoft : public Pair {
   void *extract(char *, int &);
 
  protected:
-  double PI;
   double cut_global;
   double **prefactor;
   double **cut;
diff --git a/src/set.h b/src/set.h
index 5e1fd4e8c7..00e0c77de1 100644
--- a/src/set.h
+++ b/src/set.h
@@ -35,7 +35,6 @@ class Set : protected Pointers {
   int style,ivalue,newtype,count;
   int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag;
   double dvalue,xvalue,yvalue,zvalue,wvalue,fraction;
-  double PI;
 
   void selection(int);
   void set(int);
diff --git a/src/thermo.cpp b/src/thermo.cpp
index 5bd877191a..b52a5f3c9c 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -36,10 +36,12 @@
 #include "kspace.h"
 #include "output.h"
 #include "timer.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 // customize a new keyword by adding to this list:
 
@@ -147,8 +149,6 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   format_float_user = NULL;
   format_int_user = NULL;
   format_bigint_user = NULL;
-
-  PI = 4.0*atan(1.0);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1894,7 +1894,7 @@ void Thermo::compute_cellalpha()
     double* h = domain->h;
     double cosalpha = (h[5]*h[4]+h[1]*h[3])/
       sqrt((h[1]*h[1]+h[5]*h[5])*(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]));
-    dvalue = acos(cosalpha)*180.0/PI;
+    dvalue = acos(cosalpha)*180.0/MY_PI;
   }    
 }
 
@@ -1910,7 +1910,7 @@ void Thermo::compute_cellbeta()
 
     double* h = domain->h;
     double cosbeta = h[4]/sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]);
-    dvalue = acos(cosbeta)*180.0/PI;
+    dvalue = acos(cosbeta)*180.0/MY_PI;
   }
 }
 
@@ -1926,6 +1926,6 @@ void Thermo::compute_cellgamma()
 
     double* h = domain->h;
     double cosgamma = h[5]/sqrt(h[1]*h[1]+h[5]*h[5]);
-    dvalue = acos(cosgamma)*180.0/PI;
+    dvalue = acos(cosgamma)*180.0/MY_PI;
   }
 }
diff --git a/src/thermo.h b/src/thermo.h
index 352209cb35..eeb758b8e7 100644
--- a/src/thermo.h
+++ b/src/thermo.h
@@ -96,8 +96,6 @@ class Thermo : protected Pointers {
   char **id_variable;          // list of variable names
   int *variables;              // list of Variable indices
 
-  double PI;
-
   // private methods
 
   void allocate();