use MathConst::MY_PI more consistently and benefit from it being a constexpr

This commit is contained in:
Axel Kohlmeyer
2022-03-10 14:14:32 -05:00
parent 8e85481afa
commit 0fcf40c48e
4 changed files with 24 additions and 29 deletions

View File

@ -37,6 +37,7 @@
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
enum{NONE,RLINEAR,RSQ};
@ -104,7 +105,6 @@ void PairMultiLucy::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double pi = MathConst::MY_PI;
double A_i;
double A_j;
double fraction_i,fraction_j;
@ -198,7 +198,7 @@ void PairMultiLucy::compute(int eflag, int vflag)
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
evdwl,0.0,0.0,0.0,0.0,0.0);
@ -733,7 +733,6 @@ void PairMultiLucy::computeLocalDensity()
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double pi = MathConst::MY_PI;
double *rho = atom->rho;
// zero out density
@ -766,7 +765,7 @@ void PairMultiLucy::computeLocalDensity()
if (rsq < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
factor= (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;

View File

@ -39,6 +39,7 @@
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
enum{NONE,RLINEAR,RSQ};
@ -127,7 +128,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
double pi = MathConst::MY_PI;
double A_i, A_j;
double fraction_i,fraction_j;
int jtable;
@ -276,7 +276,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwlOld = mixWtSite1old_i*evdwl;
evdwl = mixWtSite1_i*evdwl;
@ -866,15 +866,13 @@ void PairMultiLucyRX::computeLocalDensity()
const int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
const double pi = MathConst::MY_PI;
const bool newton_pair = force->newton_pair;
const bool one_type = (atom->ntypes == 1);
// Special cut-off values for when there's only one type.
const double cutsq_type11 = cutsq[1][1];
const double rcut_type11 = sqrt(cutsq_type11);
const double factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11);
const double factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11);
double *rho = atom->rho;
@ -925,7 +923,7 @@ void PairMultiLucyRX::computeLocalDensity()
const double rcut = sqrt(cutsq[itype][jtype]);
const double tmpFactor = 1.0-sqrt(rsq)/rcut;
const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho_i += factor;
if (newton_pair || j < nlocal)
rho[j] += factor;

View File

@ -23,20 +23,24 @@
------------------------------------------------------------------------------------------- */
#include "pair_multi_lucy_rx_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kokkos.h"
#include "math_const.h"
#include "math_const.h"
#include "memory_kokkos.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include <cmath>
#include <cstring>
#include "math_const.h"
#include "atom_kokkos.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory_kokkos.h"
#include "error.h"
#include "atom_masks.h"
#include "neigh_request.h"
#include "kokkos.h"
using namespace LAMMPS_NS;
using MathConst::MY_PI;
enum{NONE,RLINEAR,RSQ};
@ -282,7 +286,6 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXCompute<NEI
double mixWtSite2old_i,mixWtSite2old_j;
double mixWtSite1_i;
double pi = MathConst::MY_PI;
double A_i, A_j;
double fraction_i,fraction_j;
int jtable;
@ -417,7 +420,7 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXCompute<NEI
evdwl = d_table_const.e(tidx,itable) + fraction_i*d_table_const.de(tidx,itable);
} else k_error_flag.template view<DeviceType>()() = 3;
evdwl *=(pi*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0;
evdwl *=(MY_PI*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0;
evdwlOld = mixWtSite1old_i*evdwl;
evdwl = mixWtSite1_i*evdwl;
@ -458,15 +461,13 @@ void PairMultiLucyRXKokkos<DeviceType>::computeLocalDensity()
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
const double pi = MathConst::MY_PI;
const bool newton_pair = force->newton_pair;
const bool one_type = (atom->ntypes == 1);
// Special cut-off values for when there's only one type.
cutsq_type11 = cutsq[1][1];
rcut_type11 = sqrt(cutsq_type11);
factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11);
factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11);
// zero out density
int m = nlocal;
@ -548,8 +549,6 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXComputeLoca
const int itype = type[i];
const int jnum = d_numneigh[i];
const double pi = MathConst::MY_PI;
for (int jj = 0; jj < jnum; jj++) {
const int j = (d_neighbors(i,jj) & NEIGHMASK);
const int jtype = type[j];
@ -574,7 +573,7 @@ void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXComputeLoca
const double rcut = sqrt(d_cutsq(itype,jtype));
const double tmpFactor = 1.0-sqrt(rsq)/rcut;
const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho_i_contrib += factor;
if ((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) && (NEWTON_PAIR || j < nlocal))
a_rho[j] += factor;

View File

@ -233,7 +233,6 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3
// 21 = theta
// see alloyparams(void) in meam_setup_done.cpp
case 21:
// const double PI = 3.141592653589793238463;
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;