diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index cc92e1719f..221778e6ee 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -35,7 +35,10 @@ #include "memory.h" #include "error.h" +#include "math_const.h" + using namespace LAMMPS_NS; +using namespace MathConst; #define MAXORDER 7 #define OFFSET 16384 @@ -58,7 +61,6 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command"); precision = atof(arg[0]); - PI = 4.0*atan(1.0); nfactors = 3; factors = new int[nfactors]; @@ -520,9 +522,9 @@ void PPPM::setup() delvolinv = delxinv*delyinv*delzinv; - double unitkx = (2.0*PI/xprd); - double unitky = (2.0*PI/yprd); - double unitkz = (2.0*PI/zprd_slab); + double unitkx = (2.0*MY_PI/xprd); + double unitky = (2.0*MY_PI/yprd); + double unitkz = (2.0*MY_PI/zprd_slab); // fkx,fky,fkz for my FFT grid pts @@ -581,11 +583,11 @@ void PPPM::setup() double sum1,dot1,dot2; double numerator,denominator; - int nbx = static_cast ((g_ewald*xprd/(PI*nx_pppm)) * + int nbx = static_cast ((g_ewald*xprd/(MY_PI*nx_pppm)) * pow(-log(EPS_HOC),0.25)); - int nby = static_cast ((g_ewald*yprd/(PI*ny_pppm)) * + int nby = static_cast ((g_ewald*yprd/(MY_PI*ny_pppm)) * pow(-log(EPS_HOC),0.25)); - int nbz = static_cast ((g_ewald*zprd_slab/(PI*nz_pppm)) * + int nbz = static_cast ((g_ewald*zprd_slab/(MY_PI*nz_pppm)) * pow(-log(EPS_HOC),0.25)); double form = 1.0; @@ -708,7 +710,7 @@ void PPPM::compute(int eflag, int vflag) energy *= 0.5*volume; energy -= g_ewald*qsqsum/1.772453851 + - 0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume); + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); energy *= qqrd2e*scale; } @@ -1027,7 +1029,7 @@ double PPPM::rms(double h, double prd, bigint natoms, for (int m = 0; m < order; m++) sum += acons[order][m] * pow(h*g_ewald,2.0*m); double value = q2 * pow(h*g_ewald,order) * - sqrt(g_ewald*prd*sqrt(2.0*PI)*sum/natoms) / (prd*prd); + sqrt(g_ewald*prd*sqrt(2.0*MY_PI)*sum/natoms) / (prd*prd); return value; } @@ -1873,13 +1875,13 @@ void PPPM::slabcorr(int eflag) // compute corrections - double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume; + double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; if (eflag) energy += qqrd2e*scale * e_slabcorr; // add on force corrections - double ffact = -4.0*PI*dipole_all/volume; + double ffact = -4.0*MY_PI*dipole_all/volume; double **f = atom->f; for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact; diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index 982b8f7c9e..e2a58e8804 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -47,7 +47,6 @@ class PPPM : public KSpace { protected: int me,nprocs; - double PI; double precision; int nfactors; int *factors; diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp index 12860f7aed..5f57715f9f 100644 --- a/src/KSPACE/pppm_cg.cpp +++ b/src/KSPACE/pppm_cg.cpp @@ -26,7 +26,10 @@ #include "memory.h" #include "pppm_cg.h" +#include "math_const.h" + using namespace LAMMPS_NS; +using namespace MathConst; #define OFFSET 16384 #define SMALLQ 0.00001 @@ -176,7 +179,7 @@ void PPPMCG::compute(int eflag, int vflag) energy *= 0.5*volume; energy -= g_ewald*qsqsum/1.772453851 + - 0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume); + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); energy *= qqrd2e*scale; } @@ -372,13 +375,13 @@ void PPPMCG::slabcorr(int eflag) // compute corrections - double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume; + double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; if (eflag) energy += qqrd2e*scale * e_slabcorr; // add on force corrections - double ffact = -4.0*PI*dipole_all/volume * qqrd2e * scale; + double ffact = -4.0*MY_PI*dipole_all/volume * qqrd2e * scale; double **f = atom->f; for (int j = 0; j < num_charged; j++) { diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index bd50ec9bd3..7771bb1982 100755 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -29,7 +29,10 @@ #include "memory.h" #include "error.h" +#include "math_const.h" + using namespace LAMMPS_NS; +using namespace MathConst; #define MAXLINE 1024 #define DELTA 4 @@ -41,10 +44,6 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp) single_enable = 0; one_coeff = 1; - PI = 4.0*atan(1.0); - PI2 = 2.0*atan(1.0); - PI4 = atan(1.0); - nelements = 0; elements = NULL; nparams = maxparam = 0; @@ -634,7 +633,7 @@ double PairTersoff::ters_fc(double r, Param *param) if (r < ters_R-ters_D) return 1.0; if (r > ters_R+ters_D) return 0.0; - return 0.5*(1.0 - sin(PI2*(r - ters_R)/ters_D)); + return 0.5*(1.0 - sin(MY_PI2*(r - ters_R)/ters_D)); } /* ---------------------------------------------------------------------- */ @@ -646,7 +645,7 @@ double PairTersoff::ters_fc_d(double r, Param *param) if (r < ters_R-ters_D) return 0.0; if (r > ters_R+ters_D) return 0.0; - return -(PI4/ters_D) * cos(PI2*(r - ters_R)/ters_D); + return -(MY_PI4/ters_D) * cos(MY_PI2*(r - ters_R)/ters_D); } /* ---------------------------------------------------------------------- */ @@ -700,27 +699,6 @@ double PairTersoff::ters_bij_d(double zeta, Param *param) /* ---------------------------------------------------------------------- */ -double PairTersoff::ters_gijk(double costheta, Param *param) -{ - double ters_c = param->c; - double ters_d = param->d; - - return param->gamma*(1.0 + pow(ters_c/ters_d,2.0) - - pow(ters_c,2.0) / (pow(ters_d,2.0) + pow(param->h - costheta,2.0))); -}; - -/* ---------------------------------------------------------------------- */ - -double PairTersoff::ters_gijk_d(double costheta, Param *param) -{ - double numerator = -2.0 * pow(param->c,2) * (param->h - costheta); - double denominator = pow(pow(param->d,2.0) + - pow(param->h - costheta,2.0),2.0); - return param->gamma*numerator/denominator; -} - -/* ---------------------------------------------------------------------- */ - void PairTersoff::ters_zetaterm_d(double prefactor, double *rij_hat, double rij, double *rik_hat, double rik, @@ -794,4 +772,4 @@ void PairTersoff::costheta_d(double *rij_hat, double rij, vec3_scale(1.0/rik,drk,drk); vec3_add(drj,drk,dri); vec3_scale(-1.0,dri,dri); -}; +} diff --git a/src/MANYBODY/pair_tersoff.h b/src/MANYBODY/pair_tersoff.h index 49d3079286..92dbe635bd 100755 --- a/src/MANYBODY/pair_tersoff.h +++ b/src/MANYBODY/pair_tersoff.h @@ -28,7 +28,7 @@ class PairTersoff : public Pair { public: PairTersoff(class LAMMPS *); virtual ~PairTersoff(); - void compute(int, int); + virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); void init_style(); @@ -49,22 +49,22 @@ class PairTersoff : public Pair { double ZBLcut,ZBLexpscale; }; - double PI,PI2,PI4; - double cutmax; // max cutoff for all elements - int nelements; // # of unique elements + Param *params; // parameter set for an I-J-K interaction char **elements; // names of unique elements int ***elem2param; // mapping from element triplets to parameters int *map; // mapping from atom types to elements + double cutmax; // max cutoff for all elements + int nelements; // # of unique elements int nparams; // # of stored parameter sets int maxparam; // max # of parameter sets - Param *params; // parameter set for an I-J-K interaction void allocate(); virtual void read_file(char *); void setup(); virtual void repulsive(Param *, double, double &, int, double &); double zeta(Param *, double, double, double *, double *); - void force_zeta(Param *, double, double, double &, double &, int, double &); + virtual void force_zeta(Param *, double, double, double &, + double &, int, double &); void attractive(Param *, double, double, double, double *, double *, double *, double *, double *); @@ -74,26 +74,52 @@ class PairTersoff : public Pair { virtual double ters_fa_d(double, Param *); double ters_bij(double, Param *); double ters_bij_d(double, Param *); - double ters_gijk(double, Param *); - double ters_gijk_d(double, Param *); + void ters_zetaterm_d(double, double *, double, double *, double, double *, double *, double *, Param *); void costheta_d(double *, double, double *, double, double *, double *, double *); - // vector functions, inline for efficiency + // inlined functions for efficiency - inline double vec3_dot(double *x, double *y) { + inline double ters_gijk(const double costheta, + const Param * const param) const { + const double ters_c = param->c * param->c; + const double ters_d = param->d * param->d; + const double hcth = param->h - costheta; + + return param->gamma*(1.0 + ters_c/ters_d - ters_c / (ters_d + hcth*hcth)); + } + + inline double ters_gijk_d(const double costheta, + const Param * const param) const { + const double ters_c = param->c * param->c; + const double ters_d = param->d * param->d; + const double hcth = param->h - costheta; + const double numerator = -2.0 * ters_c * hcth; + const double denominator = 1.0/(ters_d + hcth*hcth); + return param->gamma*numerator*denominator*denominator; + } + + inline double vec3_dot(const double x[3], const double y[3]) const { return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; } - inline void vec3_add(double *x, double *y, double *z) { + + inline void vec3_add(const double x[3], const double y[3], + double * const z) const { z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2]; } - inline void vec3_scale(double k, double *x, double *y) { + + inline void vec3_scale(const double k, const double x[3], + double y[3]) const { y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2]; } - inline void vec3_scaleadd(double k, double *x, double *y, double *z) { - z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2]; + + inline void vec3_scaleadd(const double k, const double x[3], + const double y[3], double * const z) const { + z[0] = k*x[0]+y[0]; + z[1] = k*x[1]+y[1]; + z[2] = k*x[2]+y[2]; } }; diff --git a/src/MANYBODY/pair_tersoff_zbl.cpp b/src/MANYBODY/pair_tersoff_zbl.cpp index 90ccfd9c11..a4f6ee3edd 100644 --- a/src/MANYBODY/pair_tersoff_zbl.cpp +++ b/src/MANYBODY/pair_tersoff_zbl.cpp @@ -31,7 +31,9 @@ #include "memory.h" #include "error.h" +#include "math_const.h" using namespace LAMMPS_NS; +using namespace MathConst; #define MAXLINE 1024 #define DELTA 4 @@ -224,7 +226,7 @@ void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce, double esq = pow(global_e,2.0); double a_ij = (0.8854*global_a_0) / (pow(param->Z_i,0.23) + pow(param->Z_j,0.23)); - double premult = (param->Z_i * param->Z_j * esq)/(4.0*PI*global_epsilon_0); + double premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_PI*global_epsilon_0); double r_ov_a = r/a_ij; double phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) + 0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a); @@ -232,7 +234,7 @@ void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce, 0.9423*0.5099*exp(-0.9423*r_ov_a) - 0.4029*0.2802*exp(-0.4029*r_ov_a) - 0.2016*0.02817*exp(-0.2016*r_ov_a)); - double fforce_ZBL = premult*-pow(r,-2.0)* phi + premult*pow(r,-1.0)*dphi; + double fforce_ZBL = premult*-phi/rsq + premult*dphi/r; double eng_ZBL = premult*(1.0/r)*phi; // combine two parts with smoothing by Fermi-like function diff --git a/src/USER-EWALDN/math_vector.h b/src/USER-EWALDN/math_vector.h index 90fc7ee464..b7473b5d89 100644 --- a/src/USER-EWALDN/math_vector.h +++ b/src/USER-EWALDN/math_vector.h @@ -19,6 +19,7 @@ #define LMP_MATH_VECTOR_H #include "math.h" +#include "string.h" #define VECTOR_NULL {0, 0, 0} #define SHAPE_NULL {0, 0, 0, 0, 0, 0}