ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7133 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2011-10-20 00:11:49 +00:00
parent 55b951bed0
commit f90ce18f3d
83 changed files with 293 additions and 271 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/CLASS2/angle_class2.cpp
View File

@ -24,10 +24,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -318,7 +320,7 @@ void AngleClass2::coeff(int narg, char **arg)
// convert theta0 from degrees to radians // convert theta0 from degrees to radians
for (int i = ilo; i <= ihi; i++) { 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; k2[i] = k2_one;
k3[i] = k3_one; k3[i] = k3_one;
k4[i] = k4_one; k4[i] = k4_one;

21
src/CLASS2/dihedral_class2.cpp
View File

@ -26,20 +26,19 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05 #define TOLERANCE 0.05
#define SMALL 0.0000001 #define SMALL 0.0000001
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
DihedralClass2::DihedralClass2(LAMMPS *lmp) : Dihedral(lmp) DihedralClass2::DihedralClass2(LAMMPS *lmp) : Dihedral(lmp) {}
{
PI = 4.0*atan(1.0);
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -697,8 +696,8 @@ void DihedralClass2::coeff(int narg, char **arg)
at_f1_2[i] = f1_2_one; at_f1_2[i] = f1_2_one;
at_f2_2[i] = f2_2_one; at_f2_2[i] = f2_2_one;
at_f3_2[i] = f3_2_one; at_f3_2[i] = f3_2_one;
at_theta0_1[i] = theta0_1_one/180.0 * PI; at_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
at_theta0_2[i] = theta0_2_one/180.0 * PI; at_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
setflag_at[i] = 1; setflag_at[i] = 1;
count++; count++;
} }
@ -714,8 +713,8 @@ void DihedralClass2::coeff(int narg, char **arg)
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
aat_k[i] = k_one; aat_k[i] = k_one;
aat_theta0_1[i] = theta0_1_one/180.0 * PI; aat_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
aat_theta0_2[i] = theta0_2_one/180.0 * PI; aat_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
setflag_aat[i] = 1; setflag_aat[i] = 1;
count++; count++;
} }
@ -749,11 +748,11 @@ void DihedralClass2::coeff(int narg, char **arg)
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k1[i] = k1_one; k1[i] = k1_one;
phi1[i] = phi1_one/180.0 * PI; phi1[i] = phi1_one/180.0 * MY_PI;
k2[i] = k2_one; k2[i] = k2_one;
phi2[i] = phi2_one/180.0 * PI; phi2[i] = phi2_one/180.0 * MY_PI;
k3[i] = k3_one; k3[i] = k3_one;
phi3[i] = phi3_one/180.0 * PI; phi3[i] = phi3_one/180.0 * MY_PI;
setflag_d[i] = 1; setflag_d[i] = 1;
count++; count++;
} }

1
src/CLASS2/dihedral_class2.h
View File

@ -46,7 +46,6 @@ class DihedralClass2 : public Dihedral {
double *bb13t_k,*bb13t_r10,*bb13t_r30; double *bb13t_k,*bb13t_r10,*bb13t_r30;
int *setflag_d,*setflag_mbt,*setflag_ebt; int *setflag_d,*setflag_mbt,*setflag_ebt;
int *setflag_at,*setflag_aat,*setflag_bb13t; int *setflag_at,*setflag_aat,*setflag_bb13t;
double PI;
void allocate(); void allocate();
}; };

15
src/CLASS2/improper_class2.cpp
View File

@ -26,19 +26,18 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ImproperClass2::ImproperClass2(LAMMPS *lmp) : Improper(lmp) ImproperClass2::ImproperClass2(LAMMPS *lmp) : Improper(lmp) {}
{
PI = 4.0*atan(1.0);
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -550,9 +549,9 @@ void ImproperClass2::coeff(int narg, char **arg)
aa_k1[i] = k1_one; aa_k1[i] = k1_one;
aa_k2[i] = k2_one; aa_k2[i] = k2_one;
aa_k3[i] = k3_one; aa_k3[i] = k3_one;
aa_theta0_1[i] = theta0_1_one/180.0 * PI; aa_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
aa_theta0_2[i] = theta0_2_one/180.0 * PI; aa_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
aa_theta0_3[i] = theta0_3_one/180.0 * PI; aa_theta0_3[i] = theta0_3_one/180.0 * MY_PI;
setflag_aa[i] = 1; setflag_aa[i] = 1;
count++; count++;
} }
@ -567,7 +566,7 @@ void ImproperClass2::coeff(int narg, char **arg)
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k0[i] = k0_one; k0[i] = k0_one;
chi0[i] = chi0_one/180.0 * PI; chi0[i] = chi0_one/180.0 * MY_PI;
setflag_i[i] = 1; setflag_i[i] = 1;
count++; count++;
} }

1
src/CLASS2/improper_class2.h
View File

@ -38,7 +38,6 @@ class ImproperClass2 : public Improper {
double *k0,*chi0; double *k0,*chi0;
double *aa_k1,*aa_k2,*aa_k3,*aa_theta0_1,*aa_theta0_2,*aa_theta0_3; double *aa_k1,*aa_k2,*aa_k3,*aa_theta0_1,*aa_theta0_2,*aa_theta0_3;
int *setflag_i,*setflag_aa; int *setflag_i,*setflag_aa;
double PI;
void allocate(); void allocate();
void angleangle(int, int); void angleangle(int, int);

7
src/CLASS2/pair_lj_class2.cpp
View File

@ -19,10 +19,12 @@
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); 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 sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3; double sig6 = sig3*sig3;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3; 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); 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; sig6 * (sig3 - 2.0*rc3) / rc6;
} }

7
src/CLASS2/pair_lj_class2_coul_cut.cpp
View File

@ -20,10 +20,12 @@
#include "force.h" #include "force.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); 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 sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3; double sig6 = sig3*sig3;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3; 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); 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; sig6 * (sig3 - 2.0*rc3) / rc6;
} }

7
src/CLASS2/pair_lj_class2_coul_long.cpp
View File

@ -22,10 +22,12 @@
#include "kspace.h" #include "kspace.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917 #define EWALD_F 1.12837917
#define EWALD_P 0.3275911 #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); 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 sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3; double sig6 = sig3*sig3;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3; 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); 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; sig6 * (sig3 - 2.0*rc3) / rc6;
} }

15
src/COLLOID/pair_lubricate.cpp
View File

@ -28,11 +28,13 @@
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "update.h" #include "update.h"
#include "memory.h"
#include "random_mars.h" #include "random_mars.h"
#include "math_const.h"
#include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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 P_dot_wrel_1,P_dot_wrel_2,P_dot_wrel_3;
double a_squeeze,a_shear,a_pump,a_twist; double a_squeeze,a_shear,a_pump,a_twist;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double PI = 4.0*atan(1.0);
if (eflag || vflag) ev_setup(eflag,vflag); if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0; else evflag = vflag_fdotr = 0;
@ -190,16 +191,16 @@ void PairLubricate::compute(int eflag, int vflag)
h_sep = r - 2.0*radi; h_sep = r - 2.0*radi;
if (flag1) 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) 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; log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0;
if (flag3) 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) + ((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)); (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep));
if (flag4) 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)); (h_sep/2.0/radi) * log(2.0/(2.0*h_sep));
if (h_sep >= cut_inner[itype][jtype]) { if (h_sep >= cut_inner[itype][jtype]) {
@ -231,7 +232,7 @@ void PairLubricate::compute(int eflag, int vflag)
torque[i][2] += vxmu2f * tz; torque[i][2] += vxmu2f * tz;
} else { } 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]); (2.0*radi/4.0/cut_inner[itype][jtype]);
fpair = -a_squeeze*vnnr; fpair = -a_squeeze*vnnr;
fpair *= vxmu2f; fpair *= vxmu2f;

6
src/GPU/pppm_gpu.cpp
View File

@ -32,11 +32,13 @@
#include "domain.h" #include "domain.h"
#include "fft3d_wrap.h" #include "fft3d_wrap.h"
#include "remap_wrap.h" #include "remap_wrap.h"
#include "gpu_extra.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "gpu_extra.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXORDER 7 #define MAXORDER 7
#define OFFSET 16384 #define OFFSET 16384
@ -189,7 +191,7 @@ void PPPMGPU::compute(int eflag, int vflag)
energy *= 0.5*volume; energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 + 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; energy *= qqrd2e*scale;
} }

12
src/GRANULAR/fix_pour.cpp
View File

@ -27,10 +27,12 @@
#include "region_block.h" #include "region_block.h"
#include "region_cylinder.h" #include "region_cylinder.h"
#include "random_park.h" #include "random_park.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EPSILON 0.001 #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"); if (seed <= 0) error->all(FLERR,"Illegal fix pour command");
PI = 4.0*atan(1.0);
// option defaults // option defaults
int iregion = -1; int iregion = -1;
@ -214,11 +214,11 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
double dy = yhi - ylo; double dy = yhi - ylo;
if (dy < 1.0) dy = 1.0; if (dy < 1.0) dy = 1.0;
volume = (xhi-xlo) * dy * (zhi-zlo); volume = (xhi-xlo) * dy * (zhi-zlo);
} else volume = PI*rc*rc * (zhi-zlo); } else volume = MY_PI*rc*rc * (zhi-zlo);
volume_one = 4.0/3.0 * PI * radius_hi*radius_hi*radius_hi; volume_one = 4.0/3.0 * MY_PI * radius_hi*radius_hi*radius_hi;
} else { } else {
volume = (xhi-xlo) * (yhi-ylo); 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); nper = static_cast<int> (volfrac*volume/volume_one);
@ -472,7 +472,7 @@ void FixPour::pre_exchange()
m = atom->nlocal - 1; m = atom->nlocal - 1;
atom->type[m] = ntype; atom->type[m] = ntype;
atom->radius[m] = radtmp; 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->mask[m] = 1 | groupbit;
atom->v[m][0] = vxtmp; atom->v[m][0] = vxtmp;
atom->v[m][1] = vytmp; atom->v[m][1] = vytmp;

1
src/GRANULAR/fix_pour.h
View File

@ -56,7 +56,6 @@ class FixPour : public Fix {
int me,nprocs; int me,nprocs;
int *recvcounts,*displs; int *recvcounts,*displs;
double PI;
int nfreq,nfirst,ninserted,nper; int nfreq,nfirst,ninserted,nper;
double lo_current,hi_current; double lo_current,hi_current;
class FixShearHistory *fix_history; class FixShearHistory *fix_history;

7
src/GRANULAR/fix_wall_gran.cpp
View File

@ -26,10 +26,12 @@
#include "pair.h" #include "pair.h"
#include "modify.h" #include "modify.h"
#include "respa.h" #include "respa.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{XPLANE,YPLANE,ZPLANE,ZCYLINDER}; // XYZ PLANE need to be 0,1,2 enum{XPLANE,YPLANE,ZPLANE,ZCYLINDER}; // XYZ PLANE need to be 0,1,2
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY}; enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY};
@ -159,10 +161,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
// setup oscillations // setup oscillations
if (wiggle) { if (wiggle) omega = 2.0*MY_PI / period;
double PI = 4.0 * atan(1.0);
omega = 2.0*PI / period;
}
// perform initial allocation of atom-based arrays // perform initial allocation of atom-based arrays
// register with Atom class // register with Atom class

24
src/KSPACE/ewald.cpp
View File

@ -26,10 +26,12 @@
#include "force.h" #include "force.h"
#include "pair.h" #include "pair.h"
#include "domain.h" #include "domain.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.00001 #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"); if (narg != 1) error->all(FLERR,"Illegal kspace_style ewald command");
precision = atof(arg[0]); precision = atof(arg[0]);
PI = 4.0*atan(1.0);
kmax = 0; kmax = 0;
kxvecs = kyvecs = kzvecs = NULL; kxvecs = kyvecs = kzvecs = NULL;
@ -165,17 +166,17 @@ void Ewald::setup()
double zprd_slab = zprd*slab_volfactor; double zprd_slab = zprd*slab_volfactor;
volume = xprd * yprd * zprd_slab; volume = xprd * yprd * zprd_slab;
unitk[0] = 2.0*PI/xprd; unitk[0] = 2.0*MY_PI/xprd;
unitk[1] = 2.0*PI/yprd; unitk[1] = 2.0*MY_PI/yprd;
unitk[2] = 2.0*PI/zprd_slab; unitk[2] = 2.0*MY_PI/zprd_slab;
// determine kmax // determine kmax
// function of current box size, precision, G_ewald (short-range cutoff) // function of current box size, precision, G_ewald (short-range cutoff)
int nkxmx = static_cast<int> ((g_ewald*xprd/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/PI) * sqrt(-log(precision))); int nkymx = static_cast<int> ((g_ewald*yprd/MY_PI) * sqrt(-log(precision)));
int nkzmx = 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; int kmax_old = kmax;
kmax = MAX(nkxmx,nkymx); kmax = MAX(nkxmx,nkymx);
@ -281,9 +282,8 @@ void Ewald::compute(int eflag, int vflag)
for (k = 0; k < kcount; k++) for (k = 0; k < kcount; k++)
energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
sfacim_all[k]*sfacim_all[k]); sfacim_all[k]*sfacim_all[k]);
PI = 4.0*atan(1.0);
energy -= g_ewald*qsqsum/1.772453851 + 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; energy *= qqrd2e*scale;
} }
@ -495,7 +495,7 @@ void Ewald::coeffs()
double unitky = unitk[1]; double unitky = unitk[1];
double unitkz = unitk[2]; double unitkz = unitk[2];
double g_ewald_sq_inv = 1.0 / (g_ewald*g_ewald); 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; kcount = 0;
@ -817,13 +817,13 @@ void Ewald::slabcorr(int eflag)
// compute corrections // 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; if (eflag) energy += qqrd2e*scale * e_slabcorr;
// add on force corrections // 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; double **f = atom->f;
for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact; for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;

1
src/KSPACE/ewald.h
View File

@ -34,7 +34,6 @@ class Ewald : public KSpace {
double memory_usage(); double memory_usage();
private: private:
double PI;
double precision; double precision;
int kcount,kmax,kmax3d,kmax_created; int kcount,kmax,kmax3d,kmax_created;
double qqrd2e; double qqrd2e;

7
src/KSPACE/pair_born_coul_long.cpp
View File

@ -26,10 +26,12 @@
#include "kspace.h" #include "kspace.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917 #define EWALD_F 1.12837917
#define EWALD_P 0.3275911 #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); MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double PI = 4.0*atan(1.0);
double rho1 = rho[i][j]; double rho1 = rho[i][j];
double rho2 = rho1*rho1; double rho2 = rho1*rho1;
double rho3 = rho2*rho1; double rho3 = rho2*rho1;
@ -321,11 +322,11 @@ double PairBornCoulLong::init_one(int i, int j)
double rc2 = rc*rc; double rc2 = rc*rc;
double rc3 = rc2*rc; double rc3 = rc2*rc;
double rc5 = rc3*rc2; 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* (a[i][j]*exp((sigma[i][j]-rc)/rho1)*rho1*
(rc2 + 2.0*rho1*rc + 2.0*rho2) - (rc2 + 2.0*rho1*rc + 2.0*rho2) -
c[i][j]/(3.0*rc3) + d[i][j]/(5.0*rc5)); 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) * (-a[i][j]*exp((sigma[i][j]-rc)/rho1) *
(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + (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)); 2.0*c[i][j]/rc3 - 8.0*d[i][j]/(5.0*rc5));

7
src/KSPACE/pair_buck_coul_long.cpp
View File

@ -22,10 +22,12 @@
#include "kspace.h" #include "kspace.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917 #define EWALD_F 1.12837917
#define EWALD_P 0.3275911 #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); MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double PI = 4.0*atan(1.0);
double rho1 = rho[i][j]; double rho1 = rho[i][j];
double rho2 = rho1*rho1; double rho2 = rho1*rho1;
double rho3 = rho2*rho1; double rho3 = rho2*rho1;
double rc = cut_lj[i][j]; double rc = cut_lj[i][j];
double rc2 = rc*rc; double rc2 = rc*rc;
double rc3 = rc2*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) - (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) -
c[i][j]/(3.0*rc3)); 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)* (-a[i][j]*exp(-rc/rho1)*
(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3); (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3);
} }

7
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

@ -30,10 +30,12 @@
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917 #define EWALD_F 1.12837917
#define EWALD_P 0.3275911 #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); 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 sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2; double sig6 = sig2*sig2*sig2;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3; double rc6 = rc3*rc3;
double rc9 = rc3*rc6; 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); 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); sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
} }

8
src/MANYBODY/pair_airebo.cpp
View File

@ -29,10 +29,12 @@
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024 #define MAXLINE 1024
#define TOL 1.0e-9 #define TOL 1.0e-9
@ -52,8 +54,6 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp)
maxpage = 0; maxpage = 0;
pages = NULL; pages = NULL;
nC = nH = 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; dX = 0.0;
} }
else { else {
cutoff = 0.5 * (1.0+cos(PI*t)); cutoff = 0.5 * (1.0+cos(MY_PI*t));
dX = (-0.5*PI*sin(PI*t)) / (Xmax-Xmin); dX = (-MY_PI2*sin(MY_PI*t)) / (Xmax-Xmin);
} }
return cutoff; return cutoff;
} }

1
src/MANYBODY/pair_airebo.h
View File

@ -46,7 +46,6 @@ class PairAIREBO : public Pair {
int npage; // current page in page list int npage; // current page in page list
int *map; // 0 (C), 1 (H), or -1 (NULL) for each type int *map; // 0 (C), 1 (H), or -1 (NULL) for each type
double PI;
double cutlj; // user-specified LJ cutoff double cutlj; // user-specified LJ cutoff
double cutljrebosq; // cut for when to compute double cutljrebosq; // cut for when to compute
// REBO neighs of ghost atoms // REBO neighs of ghost atoms

33
src/MANYBODY/pair_comb.cpp
View File

@ -29,12 +29,14 @@
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include "group.h" #include "group.h"
#include "update.h" #include "update.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024 #define MAXLINE 1024
#define DELTA 4 #define DELTA 4
@ -47,11 +49,6 @@ PairComb::PairComb(LAMMPS *lmp) : Pair(lmp)
single_enable = 0; single_enable = 0;
one_coeff = 1; one_coeff = 1;
PI = 4.0*atan(1.0);
PI2 = 2.0*atan(1.0);
PI4 = atan(1.0);
PIsq = sqrt(PI);
nmax = 0; nmax = 0;
NCo = NULL; NCo = NULL;
bbij = NULL; bbij = NULL;
@ -930,7 +927,7 @@ double PairComb::elp(Param *param, double rsqij, double rsqik,
double rij,rik,costheta,lp1,lp3,lp6; double rij,rik,costheta,lp1,lp3,lp6;
double rmu,rmu2,comtt,fck; double rmu,rmu2,comtt,fck;
double pplp1 = param->plp1, pplp3 = param->plp3, pplp6 = param->plp6; 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 // cos(theta) of the i-j-k
// cutoff function of rik // cutoff function of rik
@ -984,7 +981,7 @@ void PairComb::flp(Param *param, double rsqij, double rsqik,
double pplp1 = param->plp1; double pplp1 = param->plp1;
double pplp3 = param->plp3; double pplp3 = param->plp3;
double pplp6 = param->plp6; 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 // 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 1.0;
if (r > comb_R+comb_D) return 0.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;
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_R) return 0.0;
if (r > comb_D) return 1.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_R) return 0.0;
if (r > comb_D) 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_R) return 1.0;
if (r > comb_D) return 0.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_R) return 0.0;
if (r > comb_D) 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; alf = 0.20;
esucon = force->qqr2e; 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; exp(-alf*alf*rcoul*rcoul)/rcoul)*esucon/rcoul;
calc3 = (erfc(rcoul*alf)/rcoul)*esucon; 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); r = sqrt(rsq);
r3 = r * rsq; r3 = r * rsq;
alf = 0.20; alf = 0.20;
alfdpi = 2.0*alf/PIsq; alfdpi = 2.0*alf/MY_PIS;
esucon = force->qqr2e; esucon = force->qqr2e;
pot_tmp = 0.0; pot_tmp = 0.0;
pot_d = 0.0; pot_d = 0.0;

1
src/MANYBODY/pair_comb.h
View File

@ -58,7 +58,6 @@ class PairComb : public Pair {
int powermint; int powermint;
}; };
double PI,PI2,PI4,PIsq;
double cutmax; // max cutoff for all elements double cutmax; // max cutoff for all elements
int nelements; // # of unique elements int nelements; // # of unique elements
char **elements; // names of unique elements char **elements; // names of unique elements

5
src/MC/fix_gcmc.cpp
View File

@ -30,10 +30,12 @@
#include "random_park.h" #include "random_park.h"
#include "force.h" #include "force.h"
#include "pair.h" #include "pair.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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 // compute beta, lambda, sigma, and the zz factor
beta = 1.0/(force->boltz*reservoir_temperature); beta = 1.0/(force->boltz*reservoir_temperature);
double PI = 4.0*atan(1.0);
double gas_mass = atom->mass[ntype]; double gas_mass = atom->mass[ntype];
double lambda = sqrt(force->hplanck*force->hplanck/ 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)); force->boltz*reservoir_temperature));
sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e); sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e);
zz = exp(beta*chemical_potential)/(pow(lambda,3)); zz = exp(beta*chemical_potential)/(pow(lambda,3));

4
src/MOLECULE/angle_charmm.cpp
View File

@ -23,10 +23,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -209,7 +211,7 @@ void AngleCharmm::coeff(int narg, char **arg)
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; 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; k_ub[i] = k_ub_one;
r_ub[i] = r_ub_one; r_ub[i] = r_ub_one;
setflag[i] = 1; setflag[i] = 1;

4
src/MOLECULE/angle_cosine.cpp
View File

@ -19,10 +19,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -174,7 +176,7 @@ void AngleCosine::coeff(int narg, char **arg)
double AngleCosine::equilibrium_angle(int i) double AngleCosine::equilibrium_angle(int i)
{ {
return PI; return MY_PI;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

4
src/MOLECULE/angle_cosine_periodic.cpp
View File

@ -23,10 +23,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -226,7 +228,7 @@ void AngleCosinePeriodic::coeff(int narg, char **arg)
double AngleCosinePeriodic::equilibrium_angle(int i) double AngleCosinePeriodic::equilibrium_angle(int i)
{ {
return PI; return MY_PI;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

4
src/MOLECULE/angle_cosine_squared.cpp
View File

@ -23,10 +23,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -178,7 +180,7 @@ void AngleCosineSquared::coeff(int narg, char **arg)
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
theta0[i] = theta0_one/180.0 * PI; theta0[i] = theta0_one/180.0 * MY_PI;
setflag[i] = 1; setflag[i] = 1;
count++; count++;
} }

4
src/MOLECULE/angle_harmonic.cpp
View File

@ -19,10 +19,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -178,7 +180,7 @@ void AngleHarmonic::coeff(int narg, char **arg)
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
theta0[i] = theta0_one/180.0 * PI; theta0[i] = theta0_one/180.0 * MY_PI;
setflag[i] = 1; setflag[i] = 1;
count++; count++;
} }

12
src/MOLECULE/angle_table.cpp
View File

@ -24,10 +24,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{LINEAR,SPLINE}; enum{LINEAR,SPLINE};
@ -238,8 +240,8 @@ void AngleTable::coeff(int narg, char **arg)
// convert theta from degrees to radians // convert theta from degrees to radians
for (int i = 0; i < tb->ninput; i++){ for (int i = 0; i < tb->ninput; i++){
tb->afile[i] *= PI/180.0; tb->afile[i] *= MY_PI/180.0;
tb->ffile[i] *= 180.0/PI; tb->ffile[i] *= 180.0/MY_PI;
} }
// spline read-in and compute a,e,f vectors within table // 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 // delta = table spacing in angle for N-1 bins
int tlm1 = tablength-1; int tlm1 = tablength-1;
tb->delta = PI/ tlm1; tb->delta = MY_PI / tlm1;
tb->invdelta = 1.0/tb->delta; tb->invdelta = 1.0/tb->delta;
tb->deltasq6 = tb->delta*tb->delta / 6.0; tb->deltasq6 = tb->delta*tb->delta / 6.0;
@ -501,8 +503,8 @@ void AngleTable::param_extract(Table *tb, char *line)
tb->fplo = atof(word); tb->fplo = atof(word);
word = strtok(NULL," \t\n\r\f"); word = strtok(NULL," \t\n\r\f");
tb->fphi = atof(word); tb->fphi = atof(word);
tb->fplo *= (180.0/PI)*(180.0/PI); tb->fplo *= (180.0/MY_PI)*(180.0/MY_PI);
tb->fphi *= (180.0/PI)*(180.0/PI); tb->fphi *= (180.0/MY_PI)*(180.0/MY_PI);
} else if (strcmp(word,"EQ") == 0) { } else if (strcmp(word,"EQ") == 0) {
word = strtok(NULL," \t\n\r\f"); word = strtok(NULL," \t\n\r\f");
tb->theta0 = atof(word); tb->theta0 = atof(word);

13
src/MOLECULE/dihedral_charmm.cpp
View File

@ -27,10 +27,12 @@
#include "force.h" #include "force.h"
#include "pair.h" #include "pair.h"
#include "update.h" #include "update.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05 #define TOLERANCE 0.05
@ -344,14 +346,12 @@ void DihedralCharmm::coeff(int narg, char **arg)
if (weight_one < 0.0 || weight_one > 1.0) if (weight_one < 0.0 || weight_one > 1.0)
error->all(FLERR,"Incorrect weight arg for dihedral coefficients"); error->all(FLERR,"Incorrect weight arg for dihedral coefficients");
double PI = 4.0*atan(1.0);
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
shift[i] = shift_one; shift[i] = shift_one;
cos_shift[i] = cos(PI*shift_one/180.0); cos_shift[i] = cos(MY_PI*shift_one/180.0);
sin_shift[i] = sin(PI*shift_one/180.0); sin_shift[i] = sin(MY_PI*shift_one/180.0);
multiplicity[i] = multiplicity_one; multiplicity[i] = multiplicity_one;
weight[i] = weight_one; weight[i] = weight_one;
setflag[i] = 1; 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(&shift[1],atom->ndihedraltypes,MPI_INT,0,world);
MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,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++) { for (int i = 1; i <= atom->ndihedraltypes; i++) {
setflag[i] = 1; setflag[i] = 1;
cos_shift[i] = cos(PI*shift[i]/180.0); cos_shift[i] = cos(MY_PI*shift[i]/180.0);
sin_shift[i] = sin(PI*shift[i]/180.0); sin_shift[i] = sin(MY_PI*shift[i]/180.0);
} }
} }

6
src/MOLECULE/dihedral_helix.cpp
View File

@ -27,10 +27,12 @@
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "update.h" #include "update.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05 #define TOLERANCE 0.05
#define SMALL 0.001 #define SMALL 0.001
@ -192,9 +194,9 @@ void DihedralHelix::compute(int eflag, int vflag)
siinv = 1.0/si; siinv = 1.0/si;
p = aphi[type]*(1.0 - c) + bphi[type]*(1.0 + cos(3.0*phi)) + 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 + 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; if (eflag) edihedral = p;

4
src/MOLECULE/improper_harmonic.cpp
View File

@ -22,10 +22,12 @@
#include "domain.h" #include "domain.h"
#include "force.h" #include "force.h"
#include "update.h" #include "update.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05 #define TOLERANCE 0.05
#define SMALL 0.001 #define SMALL 0.001
@ -248,7 +250,7 @@ void ImproperHarmonic::coeff(int narg, char **arg)
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
chi[i] = chi_one/180.0 * PI; chi[i] = chi_one/180.0 * MY_PI;
setflag[i] = 1; setflag[i] = 1;
count++; count++;
} }

4
src/MOLECULE/improper_umbrella.cpp
View File

@ -25,10 +25,12 @@
#include "domain.h" #include "domain.h"
#include "force.h" #include "force.h"
#include "update.h" #include "update.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05 #define TOLERANCE 0.05
#define SMALL 0.001 #define SMALL 0.001
@ -269,7 +271,7 @@ void ImproperUmbrella::coeff(int narg, char **arg)
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
kw[i] = k_one; 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; if (w_one == 0) C[i] = 1.0;
else C[i] = kw[i]/(pow(sin(w0[i]),2)); else C[i] = kw[i]/(pow(sin(w0[i]),2));
setflag[i] = 1; setflag[i] = 1;

14
src/MOLECULE/pair_hbond_dreiding_lj.cpp
View File

@ -26,11 +26,13 @@
#include "neighbor.h" #include "neighbor.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "domain.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "domain.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
#define CHUNK 8 #define CHUNK 8
@ -44,8 +46,6 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp)
no_virial_fdotr_compute = 1; no_virial_fdotr_compute = 1;
PI = 4.0*atan(1.0);
nparams = maxparam = 0; nparams = maxparam = 0;
params = NULL; params = NULL;
@ -159,7 +159,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
ac = acos(c); 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); s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL; if (s < SMALL) s = SMALL;
@ -279,7 +279,7 @@ void PairHbondDreidingLJ::settings(int narg, char **arg)
ap_global = force->inumeric(arg[0]); ap_global = force->inumeric(arg[0]);
cut_inner_global = force->numeric(arg[1]); cut_inner_global = force->numeric(arg[1]);
cut_outer_global = force->numeric(arg[2]); 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) if (cut_inner_one>cut_outer_one)
error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
double cut_angle_one = cut_angle_global; 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 // grow params array if necessary
if (nparams == maxparam) { 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; if (c < -1.0) c = -1.0;
ac = acos(c); 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); s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL; if (s < SMALL) s = SMALL;

1
src/MOLECULE/pair_hbond_dreiding_lj.h
View File

@ -38,7 +38,6 @@ class PairHbondDreidingLJ : public Pair {
protected: protected:
double cut_inner_global,cut_outer_global,cut_angle_global; double cut_inner_global,cut_outer_global,cut_angle_global;
int ap_global; int ap_global;
double PI;
struct Param { struct Param {
double epsilon,sigma; double epsilon,sigma;

10
src/MOLECULE/pair_hbond_dreiding_morse.cpp
View File

@ -26,11 +26,13 @@
#include "neighbor.h" #include "neighbor.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "domain.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "domain.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
#define CHUNK 8 #define CHUNK 8
@ -128,7 +130,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
ac = acos(c); 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); s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL; if (s < SMALL) s = SMALL;
@ -242,7 +244,7 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg)
if (cut_inner_one>cut_outer_one) if (cut_inner_one>cut_outer_one)
error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
double cut_angle_one = cut_angle_global; 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 // 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; if (c < -1.0) c = -1.0;
ac = acos(c); 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); s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL; if (s < SMALL) s = SMALL;

12
src/SRD/fix_srd.cpp
View File

@ -34,10 +34,12 @@
#include "fix_wall_srd.h" #include "fix_wall_srd.h"
#include "random_mars.h" #include "random_mars.h"
#include "random_park.h" #include "random_park.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{SLIP,NOSLIP}; enum{SLIP,NOSLIP};
enum{SPHERE,ELLIPSOID,WALL}; enum{SPHERE,ELLIPSOID,WALL};
@ -2087,8 +2089,6 @@ int FixSRD::update_srd(int i, double dt, double *xscoll, double *vsnew,
void FixSRD::parameterize() void FixSRD::parameterize()
{ {
double PI = 4.0*atan(1.0);
// timesteps // timesteps
dt_big = update->dt; dt_big = update->dt;
@ -2194,20 +2194,20 @@ void FixSRD::parameterize()
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & biggroupbit) { if (mask[i] & biggroupbit) {
if (radius && radius[i] > 0.0) 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) { else if (ellipsoid && ellipsoid[i] >= 0) {
double *shape = ebonus[ellipsoid[i]].shape; 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 { } else {
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & biggroupbit) { if (mask[i] & biggroupbit) {
if (radius && radius[i] > 0.0) 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) { else if (ellipsoid && ellipsoid[i] >= 0) {
double *shape = ebonus[ellipsoid[i]].shape; double *shape = ebonus[ellipsoid[i]].shape;
volbig += PI*shape[0]*shape[1]; volbig += MY_PI*shape[0]*shape[1];
} }
} }
} }

4
src/USER-CG-CMM/angle_cg_cmm.cpp
View File

@ -24,10 +24,12 @@
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 0.001 #define SMALL 0.001
@ -323,7 +325,7 @@ void AngleCGCMM::coeff(int narg, char **arg)
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
// convert theta0 from degrees to radians // 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; epsilon[i] = epsilon_one;
sigma[i] = sigma_one; sigma[i] = sigma_one;
rcut[i] = rcut_one; rcut[i] = rcut_one;

7
src/USER-CG-CMM/pair_cmm_common.cpp
View File

@ -23,8 +23,10 @@
#include "string.h" #include "string.h"
#include "ctype.h" #include "ctype.h"
#include "math.h" #include "math.h"
#include "math_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 1.0e-6 #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); 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 sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2; double sig6 = sig2*sig2*sig2;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3; double rc6 = rc3*rc3;
double rc9 = rc3*rc6; 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); 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); sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
#endif #endif
} }

8
src/USER-CUDA/fix_gravity_cuda.cpp
View File

@ -30,12 +30,13 @@
#include "update.h" #include "update.h"
#include "domain.h" #include "domain.h"
#include "respa.h" #include "respa.h"
#include "error.h"
#include "cuda.h" #include "cuda.h"
#include "cuda_modify_flags.h" #include "cuda_modify_flags.h"
#include "math_const.h"
#include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{CHUTE,SPHERICAL,GRADIENT,VECTOR}; enum{CHUTE,SPHERICAL,GRADIENT,VECTOR};
@ -79,8 +80,7 @@ FixGravityCuda::FixGravityCuda(LAMMPS *lmp, int narg, char **arg) :
zdir = atof(arg[7]); zdir = atof(arg[7]);
} else error->all(FLERR,"Illegal fix gravity command"); } else error->all(FLERR,"Illegal fix gravity command");
double PI = 4.0*atan(1.0); degree2rad = MY_PI/180.0;
degree2rad = PI/180.0;
if (style == CHUTE || style == SPHERICAL || style == GRADIENT) { if (style == CHUTE || style == SPHERICAL || style == GRADIENT) {
if (domain->dimension == 3) { if (domain->dimension == 3) {

5
src/USER-CUDA/fix_shake_cuda.cpp
View File

@ -34,8 +34,10 @@
#include "error.h" #include "error.h"
#include "cuda.h" #include "cuda.h"
#include "cuda_modify_flags.h" #include "cuda_modify_flags.h"
#include "math_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define BIG 1.0e20 #define BIG 1.0e20
#define MASSDELTA 0.1 #define MASSDELTA 0.1
@ -53,7 +55,6 @@ FixShakeCuda::FixShakeCuda(LAMMPS *lmp, int narg, char **arg) :
MPI_Comm_rank(world,&me); MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs); MPI_Comm_size(world,&nprocs);
neighbor_step=true; neighbor_step=true;
PI = 4.0*atan(1.0);
virial_flag = 1; virial_flag = 1;
create_attribute = 1; create_attribute = 1;
@ -2183,7 +2184,7 @@ void FixShakeCuda::stats()
r3 = sqrt(delx*delx + dely*dely + delz*delz); r3 = sqrt(delx*delx + dely*dely + delz*delz);
angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); 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]; m = shake_type[i][2];
a_count[m]++; a_count[m]++;
a_ave[m] += angle; a_ave[m] += angle;

1
src/USER-CUDA/fix_shake_cuda.h
View File

@ -53,7 +53,6 @@ class FixShakeCuda : public Fix {
private: private:
class Cuda *cuda; class Cuda *cuda;
int me,nprocs; int me,nprocs;
double PI;
double tolerance; // SHAKE tolerance double tolerance; // SHAKE tolerance
int max_iter; // max # of SHAKE iterations int max_iter; // max # of SHAKE iterations
int output_every; // SHAKE stat output every so often int output_every; // SHAKE stat output every so often

21
src/USER-CUDA/pppm_cuda.cpp
View File

@ -61,8 +61,10 @@
#include "cuda_wrapper_cu.h" #include "cuda_wrapper_cu.h"
#include "pppm_cuda_cu.h" #include "pppm_cuda_cu.h"
#include "cuda.h" #include "cuda.h"
#include "math_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXORDER 7 #define MAXORDER 7
#define OFFSET 4096 #define OFFSET 4096
@ -109,7 +111,6 @@ PPPMCuda::PPPMCuda(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, (narg==2?1:nar
if(narg>1) if(narg>1)
precisionmodify=arg[1][0]; precisionmodify=arg[1][0];
else precisionmodify='='; else precisionmodify='=';
PI = 4.0*atan(1.0);
nfactors = 3; nfactors = 3;
factors = new int[nfactors]; factors = new int[nfactors];
@ -648,9 +649,9 @@ void PPPMCuda::setup()
delvolinv = delxinv*delyinv*delzinv; delvolinv = delxinv*delyinv*delzinv;
double unitkx = (2.0*PI/xprd); double unitkx = (2.0*MY_PI/xprd);
double unitky = (2.0*PI/yprd); double unitky = (2.0*MY_PI/yprd);
double unitkz = (2.0*PI/zprd_slab); double unitkz = (2.0*MY_PI/zprd_slab);
// fkx,fky,fkz for my FFT grid pts // fkx,fky,fkz for my FFT grid pts
Cuda_PPPM_Setup_fkxyz_vg(nx_pppm, ny_pppm,nz_pppm,unitkx,unitky,unitkz,g_ewald); 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 sum1,dot1,dot2;
double numerator,denominator; 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)); 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)); 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)); pow(-log(EPS_HOC),0.25));
Cuda_PPPM_setup_greensfn(nx_pppm,ny_pppm,nz_pppm,unitkx,unitky,unitkz,g_ewald, Cuda_PPPM_setup_greensfn(nx_pppm,ny_pppm,nz_pppm,unitkx,unitky,unitkz,g_ewald,
nbx,nby,nbz,xprd,yprd,zprd_slab); nbx,nby,nbz,xprd,yprd,zprd_slab);
@ -967,7 +968,7 @@ else
energy *= 0.5*volume; energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 + 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; energy *= qqrd2e;
} }
@ -1728,11 +1729,11 @@ void PPPMCuda::slabcorr(int eflag)
// compute corrections // 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; 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); cuda_slabcorr_force(&cuda->shared_data,ffact);
} }

1
src/USER-EFF/atom_vec_electron.h
View File

@ -60,7 +60,6 @@ class AtomVecElectron : public AtomVec {
bigint memory_usage(); bigint memory_usage();
private: private:
double PI;
int *tag,*type,*mask,*image; int *tag,*type,*mask,*image;
double **x,**v,**f; double **x,**v,**f;
int *spin; int *spin;

6
src/USER-OMP/dihedral_helix_omp.cpp
View File

@ -25,9 +25,11 @@
#include "domain.h" #include "domain.h"
#include "force.h" #include "force.h"
#include "update.h" #include "update.h"
#include "math_const.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05 #define TOLERANCE 0.05
#define SMALL 0.001 #define SMALL 0.001
@ -211,10 +213,10 @@ void DihedralHelixOMP::eval(double **f, int nfrom, int nto, int tid)
siinv = 1.0/si; siinv = 1.0/si;
pd = -aphi[type] + 3.0*bphi[type]*sin(3.0*phi)*siinv + 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)) + 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; a = pd;

4
src/USER-OMP/pair_soft_omp.cpp
View File

@ -19,8 +19,10 @@
#include "force.h" #include "force.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define SMALL 1.0e-4 #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]) { if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq); r = sqrt(rsq);
arg = PI/cut[itype][jtype]; arg = MY_PI/cut[itype][jtype];
if (r > SMALL) fpair = factor_lj * prefactor[itype][jtype] * if (r > SMALL) fpair = factor_lj * prefactor[itype][jtype] *
sin(arg*r) * arg/r; sin(arg*r) * arg/r;
else fpair = 0.0; else fpair = 0.0;

4
src/angle.cpp
View File

@ -15,10 +15,12 @@
#include "angle.h" #include "angle.h"
#include "atom.h" #include "atom.h"
#include "force.h" #include "force.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -27,8 +29,6 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
energy = 0.0; energy = 0.0;
allocated = 0; allocated = 0;
PI = 4.0*atan(1.0);
THIRD = 1.0/3.0;
maxeatom = maxvatom = 0; maxeatom = maxvatom = 0;
eatom = NULL; eatom = NULL;

2
src/angle.h
View File

@ -40,8 +40,6 @@ class Angle : protected Pointers {
virtual double memory_usage(); virtual double memory_usage();
protected: protected:
double PI,THIRD;
int evflag; int evflag;
int eflag_either,eflag_global,eflag_atom; int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom; int vflag_either,vflag_global,vflag_atom;

6
src/atom_vec_ellipsoid.cpp
View File

@ -25,10 +25,12 @@
#include "domain.h" #include "domain.h"
#include "modify.h" #include "modify.h"
#include "fix.h" #include "fix.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 10000 #define DELTA 10000
#define DELTA_BONUS 10000 #define DELTA_BONUS 10000
@ -53,8 +55,6 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) :
atom->ellipsoid_flag = 1; atom->ellipsoid_flag = 1;
atom->rmass_flag = atom->angmom_flag = atom->torque_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; nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = NULL; bonus = NULL;
} }
@ -1200,7 +1200,7 @@ void AtomVecEllipsoid::data_atom_bonus(int m, char **values)
// reset ellipsoid mass // reset ellipsoid mass
// previously stored density in rmass // 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; bonus[nlocal_bonus].ilocal = m;
ellipsoid[m] = nlocal_bonus++; ellipsoid[m] = nlocal_bonus++;

1
src/atom_vec_ellipsoid.h
View File

@ -76,7 +76,6 @@ class AtomVecEllipsoid : public AtomVec {
void set_shape(int, double, double, double); void set_shape(int, double, double, double);
private: private:
double PI;
int *tag,*type,*mask,*image; int *tag,*type,*mask,*image;
double **x,**v,**f; double **x,**v,**f;
double *rmass; double *rmass;

10
src/atom_vec_sphere.cpp
View File

@ -23,10 +23,12 @@
#include "force.h" #include "force.h"
#include "fix.h" #include "fix.h"
#include "fix_adapt.h" #include "fix_adapt.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 10000 #define DELTA 10000
@ -50,8 +52,6 @@ AtomVecSphere::AtomVecSphere(LAMMPS *lmp, int narg, char **arg) :
atom->sphere_flag = 1; atom->sphere_flag = 1;
atom->radius_flag = atom->rmass_flag = atom->omega_flag = atom->radius_flag = atom->rmass_flag = atom->omega_flag =
atom->torque_flag = 1; 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; v[nlocal][2] = 0.0;
radius[nlocal] = 0.5; 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][0] = 0.0;
omega[nlocal][1] = 0.0; omega[nlocal][1] = 0.0;
omega[nlocal][2] = 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; if (radius[nlocal] == 0.0) rmass[nlocal] = density;
else else
rmass[nlocal] = 4.0*PI/3.0 * rmass[nlocal] = 4.0*MY_PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal] * density; radius[nlocal]*radius[nlocal]*radius[nlocal] * density;
x[nlocal][0] = coord[0]; 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; if (radius[nlocal] == 0.0) rmass[nlocal] = density;
else else
rmass[nlocal] = 4.0*PI/3.0 * rmass[nlocal] = 4.0*MY_PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal] * density; radius[nlocal]*radius[nlocal]*radius[nlocal] * density;
return 2; return 2;

1
src/atom_vec_sphere.h
View File

@ -61,7 +61,6 @@ class AtomVecSphere : public AtomVec {
bigint memory_usage(); bigint memory_usage();
private: private:
double PI;
int *tag,*type,*mask,*image; int *tag,*type,*mask,*image;
double **x,**v,**f; double **x,**v,**f;
double *radius,*density,*rmass; double *radius,*density,*rmass;

5
src/compute_angle_local.cpp
View File

@ -20,10 +20,12 @@
#include "domain.h" #include "domain.h"
#include "force.h" #include "force.h"
#include "angle.h" #include "angle.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 10000 #define DELTA 10000
@ -133,7 +135,6 @@ int ComputeAngleLocal::compute_angles(int flag)
} }
Angle *angle = force->angle; Angle *angle = force->angle;
double PI = 4.0*atan(1.0);
m = n = 0; m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) { for (atom2 = 0; atom2 < nlocal; atom2++) {
@ -170,7 +171,7 @@ int ComputeAngleLocal::compute_angles(int flag)
c /= r1*r2; c /= r1*r2;
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
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) { if (eflag >= 0) {

6
src/compute_dihedral_local.cpp
View File

@ -20,10 +20,12 @@
#include "domain.h" #include "domain.h"
#include "force.h" #include "force.h"
#include "dihedral.h" #include "dihedral.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 10000 #define DELTA 10000
#define SMALL 0.001 #define SMALL 0.001
@ -128,8 +130,6 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
} }
} }
double PI = 4.0*atan(1.0);
m = n = 0; m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) { for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue; 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;
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; n += nvalues;
} }

6
src/compute_improper_local.cpp
View File

@ -20,10 +20,12 @@
#include "domain.h" #include "domain.h"
#include "force.h" #include "force.h"
#include "improper.h" #include "improper.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 10000 #define DELTA 10000
@ -129,8 +131,6 @@ int ComputeImproperLocal::compute_impropers(int flag)
} }
} }
double PI = 4.0*atan(1.0);
m = n = 0; m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) { for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue; 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;
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; n += nvalues;
} }

7
src/compute_rdf.cpp
View File

@ -28,10 +28,12 @@
#include "neigh_request.h" #include "neigh_request.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "group.h" #include "group.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -267,10 +269,9 @@ void ComputeRDF::compute_array()
// assuming J atoms are at uniform density // assuming J atoms are at uniform density
double constant,nideal,gr,ncoord,rlower,rupper; double constant,nideal,gr,ncoord,rlower,rupper;
double PI = 4.0*atan(1.0);
if (domain->dimension == 3) { 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++) { for (m = 0; m < npairs; m++) {
ncoord = 0.0; ncoord = 0.0;
@ -289,7 +290,7 @@ void ComputeRDF::compute_array()
} }
} else { } else {
constant = PI / (domain->xprd*domain->yprd); constant = MY_PI / (domain->xprd*domain->yprd);
for (m = 0; m < npairs; m++) { for (m = 0; m < npairs; m++) {
ncoord = 0.0; ncoord = 0.0;

1
src/dihedral.cpp
View File

@ -32,7 +32,6 @@ Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp)
energy = 0.0; energy = 0.0;
allocated = 0; allocated = 0;
PI = 4.0*atan(1.0);
maxeatom = maxvatom = 0; maxeatom = maxvatom = 0;
eatom = NULL; eatom = NULL;

2
src/dihedral.h
View File

@ -41,8 +41,6 @@ class Dihedral : protected Pointers {
virtual double memory_usage(); virtual double memory_usage();
protected: protected:
double PI;
int evflag; int evflag;
int eflag_either,eflag_global,eflag_atom; int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom; int vflag_either,vflag_global,vflag_atom;

32
src/dump_image.cpp
View File

@ -30,6 +30,7 @@
#include "input.h" #include "input.h"
#include "variable.h" #include "variable.h"
#include "random_mars.h" #include "random_mars.h"
#include "math_const.h"
#include "error.h" #include "error.h"
#include "memory.h" #include "memory.h"
@ -38,6 +39,7 @@
#endif #endif
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define NCOLORS 140 #define NCOLORS 140
#define NELEMENTS 109 #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"); if (binary || multiproc) error->all(FLERR,"Invalid dump image filename");
PI = 4.0*atan(1.0);
// set filetype based on filename suffix // set filetype based on filename suffix
int n = strlen(filename); int n = strlen(filename);
@ -95,8 +95,8 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
bdiamvalue = 0.5; bdiamvalue = 0.5;
} }
width = height = 512; width = height = 512;
theta = 60.0 * PI/180.0; theta = 60.0 * MY_PI/180.0;
phi = 30.0 * PI/180.0; phi = 30.0 * MY_PI/180.0;
thetastr = phistr = NULL; thetastr = phistr = NULL;
cflag = STATIC; cflag = STATIC;
cx = cy = cz = 0.5; cx = cy = cz = 0.5;
@ -171,7 +171,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
theta = atof(arg[iarg+1]); theta = atof(arg[iarg+1]);
if (theta < 0.0 || theta > 180.0) if (theta < 0.0 || theta > 180.0)
error->all(FLERR,"Invalid dump image theta value"); 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]) { if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
int n = strlen(&arg[iarg+2][2]) + 1; 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]); strcpy(phistr,&arg[iarg+2][2]);
} else { } else {
phi = atof(arg[iarg+2]); phi = atof(arg[iarg+2]);
phi *= PI/180.0; phi *= MY_PI/180.0;
} }
iarg += 3; iarg += 3;
@ -358,25 +358,25 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
// static parameters // static parameters
FOV = PI/6.0; // 30 degrees FOV = MY_PI/6.0; // 30 degrees
ambientColor[0] = 0.0; ambientColor[0] = 0.0;
ambientColor[1] = 0.0; ambientColor[1] = 0.0;
ambientColor[2] = 0.0; ambientColor[2] = 0.0;
keyLightPhi = -PI/4.0; // -45 degrees keyLightPhi = -MY_PI4; // -45 degrees
keyLightTheta = PI/6.0; // 30 degrees keyLightTheta = MY_PI/6.0; // 30 degrees
keyLightColor[0] = 0.9; keyLightColor[0] = 0.9;
keyLightColor[1] = 0.9; keyLightColor[1] = 0.9;
keyLightColor[2] = 0.9; keyLightColor[2] = 0.9;
fillLightPhi = PI/6.0; // 30 degrees fillLightPhi = MY_PI/6.0; // 30 degrees
fillLightTheta = 0; fillLightTheta = 0;
fillLightColor[0] = 0.45; fillLightColor[0] = 0.45;
fillLightColor[1] = 0.45; fillLightColor[1] = 0.45;
fillLightColor[2] = 0.45; fillLightColor[2] = 0.45;
backLightPhi = PI; // 180 degrees backLightPhi = MY_PI; // 180 degrees
backLightTheta = PI/12.0; // 15 degrees backLightTheta = MY_PI/12.0; // 15 degrees
backLightColor[0] = 0.9; backLightColor[0] = 0.9;
backLightColor[1] = 0.9; backLightColor[1] = 0.9;
backLightColor[2] = 0.9; backLightColor[2] = 0.9;
@ -676,11 +676,11 @@ void DumpImage::view_params()
theta = input->variable->compute_equal(thetavar); theta = input->variable->compute_equal(thetavar);
if (theta < 0.0 || theta > 180.0) if (theta < 0.0 || theta > 180.0)
error->all(FLERR,"Invalid dump image theta value"); error->all(FLERR,"Invalid dump image theta value");
theta *= PI/180.0; theta *= MY_PI/180.0;
} }
if (phistr) { if (phistr) {
phi = input->variable->compute_equal(phivar); phi = input->variable->compute_equal(phivar);
phi *= PI/180.0; phi *= MY_PI/180.0;
} }
camDir[0] = sin(theta)*cos(phi); camDir[0] = sin(theta)*cos(phi);
@ -764,7 +764,7 @@ void DumpImage::view_params()
if (ssao) { if (ssao) {
SSAORadius = maxdel * 0.05 * ssaoint; SSAORadius = maxdel * 0.05 * ssaoint;
SSAOSamples = static_cast<int> (8.0 + 32.0*ssaoint); SSAOSamples = static_cast<int> (8.0 + 32.0*ssaoint);
SSAOJitter = PI / 12; SSAOJitter = MY_PI / 12;
ambientColor[0] = 0.5; ambientColor[0] = 0.5;
ambientColor[1] = 0.5; ambientColor[1] = 0.5;
ambientColor[2] = 0.5; ambientColor[2] = 0.5;
@ -1358,7 +1358,7 @@ void DumpImage::compute_SSAO()
{ {
// used for rasterizing the spheres // used for rasterizing the spheres
double delTheta = 2.0*PI / SSAOSamples; double delTheta = 2.0*MY_PI / SSAOSamples;
// typical neighborhood value for shading // typical neighborhood value for shading

2
src/dump_image.h
View File

@ -69,8 +69,6 @@ class DumpImage : public DumpCustom {
// constant view params // constant view params
double PI;
double FOV; double FOV;
double ambientColor[3]; double ambientColor[3];

9
src/fix_adapt.cpp
View File

@ -24,10 +24,12 @@
#include "kspace.h" #include "kspace.h"
#include "input.h" #include "input.h"
#include "variable.h" #include "variable.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{PAIR,KSPACE,ATOM}; enum{PAIR,KSPACE,ATOM};
enum{DIAMETER}; enum{DIAMETER};
@ -309,7 +311,6 @@ void FixAdapt::change_settings()
if (ad->aparam == DIAMETER) { if (ad->aparam == DIAMETER) {
int mflag = 0; int mflag = 0;
if (atom->rmass_flag) mflag = 1; if (atom->rmass_flag) mflag = 1;
double PI = 4.0*atan(1.0);
double density; double density;
double *radius = atom->radius; double *radius = atom->radius;
@ -324,9 +325,11 @@ void FixAdapt::change_settings()
} else { } else {
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { 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; 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;
} }
} }
} }

4
src/fix_deform.cpp
View File

@ -27,10 +27,12 @@
#include "lattice.h" #include "lattice.h"
#include "force.h" #include "force.h"
#include "modify.h" #include "modify.h"
#include "math_const.h"
#include "kspace.h" #include "kspace.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE}; enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE};
enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; 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); if (force_reneighbor) irregular = new Irregular(lmp);
else irregular = NULL; else irregular = NULL;
TWOPI = 8.0*atan(1.0); TWOPI = 2.0*MY_PI;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

5
src/fix_gravity.cpp
View File

@ -20,9 +20,11 @@
#include "update.h" #include "update.h"
#include "domain.h" #include "domain.h"
#include "respa.h" #include "respa.h"
#include "math_const.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{CHUTE,SPHERICAL,GRADIENT,VECTOR}; enum{CHUTE,SPHERICAL,GRADIENT,VECTOR};
@ -65,8 +67,7 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
zdir = atof(arg[7]); zdir = atof(arg[7]);
} else error->all(FLERR,"Illegal fix gravity command"); } else error->all(FLERR,"Illegal fix gravity command");
double PI = 4.0*atan(1.0); degree2rad = MY_PI/180.0;
degree2rad = PI/180.0;
if (style == CHUTE || style == SPHERICAL || style == GRADIENT) { if (style == CHUTE || style == SPHERICAL || style == GRADIENT) {
if (domain->dimension == 3) { if (domain->dimension == 3) {

7
src/fix_move.cpp
View File

@ -26,10 +26,12 @@
#include "respa.h" #include "respa.h"
#include "input.h" #include "input.h"
#include "variable.h" #include "variable.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{LINEAR,WIGGLE,ROTATE,VARIABLE}; enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
enum{EQUAL,ATOM}; enum{EQUAL,ATOM};
@ -216,10 +218,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
// set omega_rotate from period // set omega_rotate from period
if (mstyle == WIGGLE || mstyle == ROTATE) { if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = 2.0*MY_PI / period;
double PI = 4.0 * atan(1.0);
omega_rotate = 2.0*PI / period;
}
// runit = unit vector along rotation axis // runit = unit vector along rotation axis

9
src/fix_orient_fcc.cpp
View File

@ -28,10 +28,12 @@
#include "neigh_request.h" #include "neigh_request.h"
#include "comm.h" #include "comm.h"
#include "output.h" #include "output.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define BIG 1000000000 #define BIG 1000000000
@ -77,7 +79,6 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
// initializations // initializations
PI = 4.0*atan(1.0);
half_fcc_nn = 6; half_fcc_nn = 6;
use_xismooth = false; use_xismooth = false;
double xicutoff = 1.57; double xicutoff = 1.57;
@ -344,10 +345,10 @@ void FixOrientFCC::post_force(int vflag)
edelta = Vxi; edelta = Vxi;
order[i][1] = 1.0; order[i][1] = 1.0;
} else { } else {
omega = (0.5*PI)*(xi_total-xi0) / (xi1-xi0); omega = MY_PI2*(xi_total-xi0) / (xi1-xi0);
nbr[i].duxi = PI*Vxi*sin(2.0*omega) / (2.0*(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; 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; added_energy += edelta;
} }

1
src/fix_orient_fcc.h
View File

@ -58,7 +58,6 @@ class FixOrientFCC : public Fix {
private: private:
int me; int me;
double PI;
int nlevels_respa; int nlevels_respa;
int direction_of_motion; // 1 = center shrinks, 0 = center grows int direction_of_motion; // 1 = center shrinks, 0 = center grows

5
src/fix_restrain.cpp
View File

@ -26,10 +26,12 @@
#include "comm.h" #include "comm.h"
#include "respa.h" #include "respa.h"
#include "input.h" #include "input.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
enum{DIHEDRAL}; enum{DIHEDRAL};
@ -83,13 +85,12 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
// special treatment for dihedral restraints // special treatment for dihedral restraints
if (rstyle == DIHEDRAL) { if (rstyle == DIHEDRAL) {
double PI = 4.0*atan(1.0);
cos_shift = (double *) cos_shift = (double *)
memory->smalloc(n_bonds * sizeof(double), "restrain:cos_shift"); memory->smalloc(n_bonds * sizeof(double), "restrain:cos_shift");
sin_shift = (double *) sin_shift = (double *)
memory->smalloc(n_bonds * sizeof(double), "restrain:sin_shift"); memory->smalloc(n_bonds * sizeof(double), "restrain:sin_shift");
for (ibond = 0; ibond < n_bonds; ibond++) { 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); cos_shift[ibond] = cos(my_arg);
sin_shift[ibond] = sin(my_arg); sin_shift[ibond] = sin(my_arg);
} }

6
src/fix_shake.cpp
View File

@ -30,10 +30,12 @@
#include "comm.h" #include "comm.h"
#include "group.h" #include "group.h"
#include "fix_respa.h" #include "fix_respa.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define BIG 1.0e20 #define BIG 1.0e20
#define MASSDELTA 0.1 #define MASSDELTA 0.1
@ -46,8 +48,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
MPI_Comm_rank(world,&me); MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs); MPI_Comm_size(world,&nprocs);
PI = 4.0*atan(1.0);
virial_flag = 1; virial_flag = 1;
create_attribute = 1; create_attribute = 1;
@ -2114,7 +2114,7 @@ void FixShake::stats()
r3 = sqrt(delx*delx + dely*dely + delz*delz); r3 = sqrt(delx*delx + dely*dely + delz*delz);
angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); 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]; m = shake_type[i][2];
a_count[m]++; a_count[m]++;
a_ave[m] += angle; a_ave[m] += angle;

1
src/fix_shake.h
View File

@ -49,7 +49,6 @@ class FixShake : public Fix {
private: private:
int me,nprocs; int me,nprocs;
double PI;
double tolerance; // SHAKE tolerance double tolerance; // SHAKE tolerance
int max_iter; // max # of SHAKE iterations int max_iter; // max # of SHAKE iterations
int output_every; // SHAKE stat output every so often int output_every; // SHAKE stat output every so often

1
src/improper.cpp
View File

@ -27,7 +27,6 @@ Improper::Improper(LAMMPS *lmp) : Pointers(lmp)
energy = 0.0; energy = 0.0;
allocated = 0; allocated = 0;
PI = 4.0*atan(1.0);
maxeatom = maxvatom = 0; maxeatom = maxvatom = 0;
eatom = NULL; eatom = NULL;

2
src/improper.h
View File

@ -38,8 +38,6 @@ class Improper : protected Pointers {
virtual double memory_usage(); virtual double memory_usage();
protected: protected:
double PI;
int evflag; int evflag;
int eflag_either,eflag_global,eflag_atom; int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom; int vflag_either,vflag_global,vflag_atom;

7
src/pair_born.cpp
View File

@ -24,10 +24,12 @@
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double PI = 4.0*atan(1.0);
double rho1 = rho[i][j]; double rho1 = rho[i][j];
double rho2 = rho1*rho1; double rho2 = rho1*rho1;
double rho3 = rho2*rho1; double rho3 = rho2*rho1;
@ -277,11 +278,11 @@ double PairBorn::init_one(int i, int j)
double rc2 = rc*rc; double rc2 = rc*rc;
double rc3 = rc2*rc; double rc3 = rc2*rc;
double rc5 = rc3*rc2; 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* (a[i][j]*exp((sigma[i][j]-rc)/rho1)*rho1*
(rc2 + 2.0*rho1*rc + 2.0*rho2) - (rc2 + 2.0*rho1*rc + 2.0*rho2) -
c[i][j]/(3.0*rc3) + d[i][j]/(5.0*rc5)); 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) * (-a[i][j]*exp((sigma[i][j]-rc)/rho1) *
(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + (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)); 2.0*c[i][j]/rc3 - 8.0*d[i][j]/(5.0*rc5));

7
src/pair_buck.cpp
View File

@ -20,10 +20,12 @@
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double PI = 4.0*atan(1.0);
double rho1 = rho[i][j]; double rho1 = rho[i][j];
double rho2 = rho1*rho1; double rho2 = rho1*rho1;
double rho3 = rho2*rho1; double rho3 = rho2*rho1;
double rc = cut[i][j]; double rc = cut[i][j];
double rc2 = rc*rc; double rc2 = rc*rc;
double rc3 = rc2*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) - (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) -
c[i][j]/(3.0*rc3)); 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)* (-a[i][j]*exp(-rc/rho1)*
(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3); (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3);
} }

7
src/pair_buck_coul_cut.cpp
View File

@ -24,10 +24,12 @@
#include "force.h" #include "force.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double PI = 4.0*atan(1.0);
double rho1 = rho[i][j]; double rho1 = rho[i][j];
double rho2 = rho1*rho1; double rho2 = rho1*rho1;
double rho3 = rho2*rho1; double rho3 = rho2*rho1;
double rc = cut_lj[i][j]; double rc = cut_lj[i][j];
double rc2 = rc*rc; double rc2 = rc*rc;
double rc3 = rc2*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) - (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) -
c[i][j]/(3.0*rc3)); 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)* (-a[i][j]*exp(-rc/rho1)*
(rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3); (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3);
} }

7
src/pair_lj96_cut.cpp
View File

@ -29,10 +29,12 @@
#include "update.h" #include "update.h"
#include "integrate.h" #include "integrate.h"
#include "respa.h" #include "respa.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); 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 sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3; double sig6 = sig3*sig3;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3; 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); 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); sig6 * (3.0*sig3 - 4.0*rc3) / (6.0*rc6);
} }

7
src/pair_lj_cut.cpp
View File

@ -29,10 +29,12 @@
#include "update.h" #include "update.h"
#include "integrate.h" #include "integrate.h"
#include "respa.h" #include "respa.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); 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 sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2; double sig6 = sig2*sig2*sig2;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3; double rc6 = rc3*rc3;
double rc9 = rc3*rc6; 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); 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); sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
} }

7
src/pair_lj_cut_coul_cut.cpp
View File

@ -21,10 +21,12 @@
#include "force.h" #include "force.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); 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 sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2; double sig6 = sig2*sig2*sig2;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3; double rc6 = rc3*rc3;
double rc9 = rc3*rc6; 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); 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); sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
} }

7
src/pair_lj_expand.cpp
View File

@ -19,10 +19,12 @@
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; 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); 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 sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2; double sig6 = sig2*sig2*sig2;
double shiftcut = shift[i][j] - cut[i][j]; 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 rc12 = rc11*shiftcut;
double shift2 = shift[i][j]*shift[i][j]; double shift2 = shift[i][j]*shift[i][j];
double shift3 = shift2*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) - sig6*((-1.0/(9.0*rc9) + shift[i][j]/(5.0*rc10) -
shift2/(11.0*rc11))*sig6 + shift2/(11.0*rc11))*sig6 +
1.0/(3.0*rc3) - shift[i][j]/(2.0*rc4) + shift2/(5.0*rc5)); 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) - sig6* ((-4.0/(3.0*rc9) + 18.0*shift[i][j]/(5.0*rc10) -
36.0*shift2/(11.0*rc11) + shift3/rc12)*sig6 + 36.0*shift2/(11.0*rc11) + shift3/rc12)*sig6 +
2.0/rc3 - 9.0*shift[i][j]/(2.0*rc4) + 2.0/rc3 - 9.0*shift[i][j]/(2.0*rc4) +

15
src/pair_soft.cpp
View File

@ -21,17 +21,16 @@
#include "force.h" #include "force.h"
#include "update.h" #include "update.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
PairSoft::PairSoft(LAMMPS *lmp) : Pair(lmp) PairSoft::PairSoft(LAMMPS *lmp) : Pair(lmp) {}
{
PI = 4.0*atan(1.0);
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -95,9 +94,9 @@ void PairSoft::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq); 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] * 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; else fpair = 0.0;
f[i][0] += delx*fpair; 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; double r,arg,philj;
r = sqrt(rsq); r = sqrt(rsq);
arg = PI*r/cut[itype][jtype]; arg = MY_PI*r/cut[itype][jtype];
fforce = factor_lj * prefactor[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)); philj = prefactor[itype][jtype] * (1.0+cos(arg));
return factor_lj*philj; return factor_lj*philj;

1
src/pair_soft.h
View File

@ -43,7 +43,6 @@ class PairSoft : public Pair {
void *extract(char *, int &); void *extract(char *, int &);
protected: protected:
double PI;
double cut_global; double cut_global;
double **prefactor; double **prefactor;
double **cut; double **cut;

1
src/set.h
View File

@ -35,7 +35,6 @@ class Set : protected Pointers {
int style,ivalue,newtype,count; int style,ivalue,newtype,count;
int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag; int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag;
double dvalue,xvalue,yvalue,zvalue,wvalue,fraction; double dvalue,xvalue,yvalue,zvalue,wvalue,fraction;
double PI;
void selection(int); void selection(int);
void set(int); void set(int);

10
src/thermo.cpp
View File

@ -36,10 +36,12 @@
#include "kspace.h" #include "kspace.h"
#include "output.h" #include "output.h"
#include "timer.h" #include "timer.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
// customize a new keyword by adding to this list: // 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_float_user = NULL;
format_int_user = NULL; format_int_user = NULL;
format_bigint_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* h = domain->h;
double cosalpha = (h[5]*h[4]+h[1]*h[3])/ 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])); 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* h = domain->h;
double cosbeta = h[4]/sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]); 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* h = domain->h;
double cosgamma = h[5]/sqrt(h[1]*h[1]+h[5]*h[5]); 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;
} }
} }

2
src/thermo.h
View File

@ -96,8 +96,6 @@ class Thermo : protected Pointers {
char **id_variable; // list of variable names char **id_variable; // list of variable names
int *variables; // list of Variable indices int *variables; // list of Variable indices
double PI;
// private methods // private methods
void allocate(); void allocate();