diff --git a/src/FLD/pair_brownian.cpp b/src/FLD/pair_brownian.cpp index 6ac2ffea3e..586c83cd0a 100755 --- a/src/FLD/pair_brownian.cpp +++ b/src/FLD/pair_brownian.cpp @@ -129,11 +129,11 @@ void PairBrownian::compute(int eflag, int vflag) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); + RT0 = 8*MY_PI*mu*pow(rad,3.0); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -208,7 +208,7 @@ void PairBrownian::compute(int eflag, int vflag) if (flaglog) { a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); @@ -541,7 +541,7 @@ void PairBrownian::init_style() // vol_P = volume of particles, assuming mono-dispersity // vol_f = volume fraction - vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3); + vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3.0); double vol_f = vol_P/vol_T; @@ -550,10 +550,10 @@ void PairBrownian::init_style() if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); // not actually needed + RT0 = 8*MY_PI*mu*pow(rad,3.0); // not actually needed } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); } } @@ -717,7 +717,7 @@ void PairBrownian::set_3_orthogonal_vectors(double p1[3], // normalize p2 - norm = sqrt(pow(p2[0],2) + pow(p2[1],2) + pow(p2[2],2)); + norm = sqrt(p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2]); p2[0] = p2[0]/norm; p2[1] = p2[1]/norm; diff --git a/src/FLD/pair_brownian_poly.cpp b/src/FLD/pair_brownian_poly.cpp index 437d2c6714..eeec6e33fc 100644 --- a/src/FLD/pair_brownian_poly.cpp +++ b/src/FLD/pair_brownian_poly.cpp @@ -115,11 +115,11 @@ void PairBrownianPoly::compute(int eflag, int vflag) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); + RT0 = 8*MY_PI*mu*pow(rad,3.0); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -151,11 +151,12 @@ void PairBrownianPoly::compute(int eflag, int vflag) f[i][1] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5); f[i][2] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5); if (flaglog) { - torque[i][0] += prethermostat*sqrt(RT0*pow(radi,3)) * + const double radi3 = radi*radi*radi; + torque[i][0] += prethermostat*sqrt(RT0*radi3) * (random->uniform()-0.5); - torque[i][1] += prethermostat*sqrt(RT0*pow(radi,3)) * + torque[i][1] += prethermostat*sqrt(RT0*radi3) * (random->uniform()-0.5); - torque[i][2] += prethermostat*sqrt(RT0*pow(radi,3)) * + torque[i][2] += prethermostat*sqrt(RT0*radi3) * (random->uniform()-0.5); } } @@ -199,20 +200,20 @@ void PairBrownianPoly::compute(int eflag, int vflag) if (flaglog) { a_sq = beta0*beta0/beta1/beta1/h_sep + - (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); - a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) + - pow(beta0,4))/21.0/pow(beta1,4)*h_sep*log(1.0/h_sep); + (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3.0)*log(1.0/h_sep); + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3.0) + + pow(beta0,4.0))/21.0/pow(beta1,4.0)*h_sep*log(1.0/h_sep); a_sq *= 6.0*MY_PI*mu*radi; - a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) * + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3.0) * log(1.0/h_sep); - a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + - 16.0*pow(beta0,4))/375.0/pow(beta1,4) * + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3.0) + + 16.0*pow(beta0,4.0))/375.0/pow(beta1,4.0) * h_sep*log(1.0/h_sep); a_sh *= 6.0*MY_PI*mu*radi; a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 * - pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); - a_pu *= 8.0*MY_PI*mu*pow(radi,3); + pow(beta0,3.0))/250.0/pow(beta1,3.0)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3.0); } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); @@ -383,7 +384,7 @@ void PairBrownianPoly::init_style() double volP = 0.0; for (int i = 0; i < nlocal; i++) - volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3); + volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3.0); MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world); double vol_f = vol_P/vol_T; diff --git a/src/FLD/pair_lubricate.cpp b/src/FLD/pair_lubricate.cpp index aca4306aaa..53f348b8b9 100755 --- a/src/FLD/pair_lubricate.cpp +++ b/src/FLD/pair_lubricate.cpp @@ -196,13 +196,13 @@ void PairLubricate::compute(int eflag, int vflag) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)* + RT0 = 8*MY_PI*mu*pow(rad,3.0); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)* (1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)* + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)* (1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -320,7 +320,7 @@ void PairLubricate::compute(int eflag, int vflag) if (flaglog) { a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); @@ -638,7 +638,7 @@ void PairLubricate::init_style() // vol_P = volume of particles, assuming monodispersity // vol_f = volume fraction - vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3); + vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3.0); double vol_f = vol_P/vol_T; if (!flagVF) vol_f = 0; @@ -647,12 +647,12 @@ void PairLubricate::init_style() if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } diff --git a/src/FLD/pair_lubricateU.cpp b/src/FLD/pair_lubricateU.cpp index c8b6f7fea2..419ac3fcbe 100644 --- a/src/FLD/pair_lubricateU.cpp +++ b/src/FLD/pair_lubricateU.cpp @@ -642,12 +642,12 @@ void PairLubricateU::compute_Fh(double **x) if (flaglog == 0) { // R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); // RT0 = 8*MY_PI*mu*pow(rad,3); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)* + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)* (1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { // R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); // RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)* + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)* (1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -880,11 +880,11 @@ void PairLubricateU::compute_RU() double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); + RT0 = 8*MY_PI*mu*pow(rad,3.0); // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -989,7 +989,7 @@ void PairLubricateU::compute_RU() if (flaglog) { a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); @@ -1157,11 +1157,11 @@ void PairLubricateU::compute_RU(double **x) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); + RT0 = 8*MY_PI*mu*pow(rad,3.0); // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -1266,7 +1266,7 @@ void PairLubricateU::compute_RU(double **x) if (flaglog) { a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); @@ -1467,7 +1467,7 @@ void PairLubricateU::compute_RE() if (flaglog) { a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } // Relative velocity at the point of closest approach due to Ef only @@ -1648,7 +1648,7 @@ void PairLubricateU::compute_RE(double **x) if (flaglog) { a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } // Relative velocity at the point of closest approach due to Ef only @@ -1921,7 +1921,7 @@ void PairLubricateU::init_style() if (atom->radius) tmp = atom->radius[0]; MPI_Allreduce(&tmp,&rad,1,MPI_DOUBLE,MPI_MAX,world); - vol_P = atom->natoms * (4.0/3.0)*MY_PI*pow(rad,3); + vol_P = atom->natoms * (4.0/3.0)*MY_PI*pow(rad,3.0); // vol_f = volume fraction @@ -1933,12 +1933,12 @@ void PairLubricateU::init_style() if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); // not actually needed - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0); // not actually needed + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } diff --git a/src/FLD/pair_lubricateU_poly.cpp b/src/FLD/pair_lubricateU_poly.cpp index b4ae44c3a1..4662910241 100644 --- a/src/FLD/pair_lubricateU_poly.cpp +++ b/src/FLD/pair_lubricateU_poly.cpp @@ -434,7 +434,7 @@ void PairLubricateUPoly::compute_Fh(double **x) // need to set delx and fy only fx = 0.0; delx = radi; - fy = vxmu2f*RS0*pow(radi,3)*gdot/2.0/radi; dely = 0.0; + fy = vxmu2f*RS0*pow(radi,3.0)*gdot/2.0/radi; dely = 0.0; fz = 0.0; delz = 0.0; if (evflag) ev_tally_xyz(i,i,nlocal,newton_pair,0.0,0.0,-fx,-fy,-fz,delx,dely,delz); @@ -717,9 +717,10 @@ void PairLubricateUPoly::compute_RU(double **x) f[i][1] += -vxmu2f*R0*radi*v[i][1]; f[i][2] += -vxmu2f*R0*radi*v[i][2]; - torque[i][0] += -vxmu2f*RT0*pow(radi,3)*wi[0]; - torque[i][1] += -vxmu2f*RT0*pow(radi,3)*wi[1]; - torque[i][2] += -vxmu2f*RT0*pow(radi,3)*wi[2]; + const double radi3 = radi*radi*radi; + torque[i][0] += -vxmu2f*RT0*radi3*wi[0]; + torque[i][1] += -vxmu2f*RT0*radi3*wi[1]; + torque[i][2] += -vxmu2f*RT0*radi3*wi[2]; if (!flagHI) continue; @@ -1251,7 +1252,7 @@ void PairLubricateUPoly::init_style() double volP = 0.0; for (int i = 0; i < nlocal; i++) - volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3); + volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3.0); MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world); double vol_f = vol_P/vol_T; diff --git a/src/FLD/pair_lubricate_poly.cpp b/src/FLD/pair_lubricate_poly.cpp index 21dd3117a6..f03fae3dfa 100644 --- a/src/FLD/pair_lubricate_poly.cpp +++ b/src/FLD/pair_lubricate_poly.cpp @@ -214,13 +214,14 @@ void PairLubricatePoly::compute(int eflag, int vflag) if (flagfld) { f[i][0] -= vxmu2f*R0*radi*v[i][0]; f[i][1] -= vxmu2f*R0*radi*v[i][1]; - f[i][2] -= vxmu2f*R0*radi*v[i][2]; - torque[i][0] -= vxmu2f*RT0*pow(radi,3)*wi[0]; - torque[i][1] -= vxmu2f*RT0*pow(radi,3)*wi[1]; - torque[i][2] -= vxmu2f*RT0*pow(radi,3)*wi[2]; + f[i][2] -= vxmu2f*R0*radi*v[i][2]; + const double radi3 = radi*radi*radi; + torque[i][0] -= vxmu2f*RT0*radi3*wi[0]; + torque[i][1] -= vxmu2f*RT0*radi3*wi[1]; + torque[i][2] -= vxmu2f*RT0*radi3*wi[2]; if (shearing && vflag_either) { - vRS0 = -vxmu2f * RS0*pow(radi,3); + vRS0 = -vxmu2f * RS0*radi3; v_tally_tensor(i,i,nlocal,newton_pair, vRS0*Ef[0][0],vRS0*Ef[1][1],vRS0*Ef[2][2], vRS0*Ef[0][1],vRS0*Ef[0][2],vRS0*Ef[1][2]); @@ -304,21 +305,21 @@ void PairLubricatePoly::compute(int eflag, int vflag) if (flaglog) { a_sq = beta0*beta0/beta1/beta1/h_sep + - (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); + (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3.0)*log(1.0/h_sep); a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0 * - pow(beta0,3)+pow(beta0,4))/21.0/pow(beta1,4) * + pow(beta0,3.0)+pow(beta0,4.0))/21.0/pow(beta1,4.0) * h_sep*log(1.0/h_sep); a_sq *= 6.0*MY_PI*mu*radi; - a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) * + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3.0) * log(1.0/h_sep); - a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + - 16.0*pow(beta0,4))/375.0/pow(beta1,4) * + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3.0) + + 16.0*pow(beta0,4.0))/375.0/pow(beta1,4.0) * h_sep*log(1.0/h_sep); a_sh *= 6.0*MY_PI*mu*radi; a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 * - pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); - a_pu *= 8.0*MY_PI*mu*pow(radi,3); + pow(beta0,3.0))/250.0/pow(beta1,3.0)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3.0); } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); // relative velocity at the point of closest approach @@ -527,7 +528,7 @@ void PairLubricatePoly::init_style() double volP = 0.0; for (int i = 0; i < nlocal; i++) - volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3); + volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3.0); MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world); double vol_f = vol_P/vol_T; diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 693a6015bf..8c70c87775 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -157,7 +157,7 @@ void PPPM::init() if (force->pair == NULL) error->all(FLERR,"KSpace style is incompatible with Pair style"); - int itmp; + int itmp=0; double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp); if (p_cutoff == NULL) error->all(FLERR,"KSpace style is incompatible with Pair style"); @@ -648,24 +648,25 @@ void PPPM::setup() numerator = form*12.5663706/sqk; denominator = gf_denom(snx2,sny2,snz2); sum1 = 0.0; + const double dorder = static_cast(order); for (nx = -nbx; nx <= nbx; nx++) { qx = unitkx*(kper+nx_pppm*nx); sx = exp(-0.25*pow(qx/g_ewald,2.0)); wx = 1.0; argx = 0.5*qx*xprd/nx_pppm; - if (argx != 0.0) wx = pow(sin(argx)/argx,order); + if (argx != 0.0) wx = pow(sin(argx)/argx,dorder); for (ny = -nby; ny <= nby; ny++) { qy = unitky*(lper+ny_pppm*ny); sy = exp(-0.25*pow(qy/g_ewald,2.0)); wy = 1.0; argy = 0.5*qy*yprd/ny_pppm; - if (argy != 0.0) wy = pow(sin(argy)/argy,order); + if (argy != 0.0) wy = pow(sin(argy)/argy,dorder); for (nz = -nbz; nz <= nbz; nz++) { qz = unitkz*(mper+nz_pppm*nz); sz = exp(-0.25*pow(qz/g_ewald,2.0)); wz = 1.0; argz = 0.5*qz*zprd_slab/nz_pppm; - if (argz != 0.0) wz = pow(sin(argz)/argz,order); + if (argz != 0.0) wz = pow(sin(argz)/argz,dorder); dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; dot2 = qx*qx+qy*qy+qz*qz; @@ -1051,9 +1052,9 @@ void PPPM::set_grid() // adjust g_ewald for new grid size - h_x = xprd/nx_pppm; - h_y = yprd/ny_pppm; - h_z = zprd_slab/nz_pppm; + h_x = xprd/static_cast(nx_pppm); + h_y = yprd/static_cast(ny_pppm); + h_z = zprd_slab/static_cast(nz_pppm); if (!gewaldflag) { double gew1,gew2,dgew,f,fmid,hmin,rtb; @@ -1157,7 +1158,7 @@ double PPPM::rms(double h, double prd, bigint natoms, double sum = 0.0; 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) * + double value = q2 * pow(h*g_ewald,(double)order) * sqrt(g_ewald*prd*sqrt(2.0*MY_PI)*sum/natoms) / (prd*prd); return value; } diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index f0e7718234..07d5d8bf90 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -714,7 +714,7 @@ void PairAIREBO::FLJ(int eflag, int vflag) drij = rij - rljmin; swidth = rljmax - rljmin; tee = drij / swidth; - tee2 = pow (tee,2); + tee2 = tee*tee; slw = 1.0 - tee2 * (3.0 - 2.0 * tee); dslw = 6.0 * tee * (1.0 - tee) / rij / swidth; } else { diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index f0e9a86ef0..004c96ce3f 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -790,7 +790,7 @@ void PairComb::setup() params[m].Qo1 = (params[m].QU1+params[m].QL1)/2.0; // (A22) params[m].dQ1 = (params[m].QU1-params[m].QL1)/2.0; // (A21) params[m].aB1 = 1.0 / - (1.0-pow(fabs(params[m].Qo1/params[m].dQ1),10)); // (A20) + (1.0-pow(fabs(params[m].Qo1/params[m].dQ1),10.0)); // (A20) params[m].bB1 = pow(fabs(params[m].aB1),0.1)/params[m].dQ1; // (A19) params[m].nD1 = log(params[m].DU1/(params[m].DU1-params[m].DL1))/ log(params[m].QU1/(params[m].QU1-params[m].QL1)); @@ -800,7 +800,7 @@ void PairComb::setup() params[m].Qo2 = (params[m].QU2+params[m].QL2)/2.0; // (A22) params[m].dQ2 = (params[m].QU2-params[m].QL2)/2.0; // (A21) params[m].aB2 = 1.0 / - (1.0-pow(fabs(params[m].Qo2/params[m].dQ2),10)); // (A20) + (1.0-pow(fabs(params[m].Qo2/params[m].dQ2),10.0)); // (A20) params[m].bB2 = pow(fabs(params[m].aB2),0.1)/params[m].dQ2; // (A19) params[m].nD2 = log(params[m].DU2/(params[m].DU2-params[m].DL2))/ log(params[m].QU2/(params[m].QU2-params[m].QL2)); @@ -1148,8 +1148,8 @@ double PairComb::self(Param *param, double qi, double selfpot) self_tmp = qi*(s1+qi*(s2+selfpot+qi*(s3+qi*(s4+qi*qi*s5)))); - if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4); - if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4); + if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4.0); + if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4.0); return self_tmp; } @@ -1167,9 +1167,9 @@ double PairComb::comb_fa(double r, Param *param, double iq, double jq) Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1); Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2); Bsi = param->bigb1 * exp(param->lam21*Di)* - (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10))); + (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10.0))); Bsj = param->bigb2 * exp(param->lam22*Dj)* - (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10))); + (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10.0))); if (Bsi > 0.0 && Bsj > 0.0) bigB = sqrt(Bsi*Bsj)*param->romigb; else bigB = 0.0; @@ -1189,9 +1189,9 @@ double PairComb::comb_fa_d(double r, Param *param, double iq, double jq) Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1); Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2); Bsi = param->bigb1 * exp(param->lam21*Di)* - (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10))); + (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10.0))); Bsj = param->bigb2 * exp(param->lam22*Dj)* - (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10))); + (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10.0))); if (Bsi > 0.0 && Bsj > 0.0) bigB = sqrt(Bsi*Bsj)*param->romigb; else bigB = 0.0; @@ -1206,7 +1206,7 @@ double PairComb::comb_bij(double zeta, Param *param) double tmp = param->beta * zeta; if (tmp > param->c1) return 1.0/sqrt(tmp); if (tmp > param->c2) - return (1.0 - pow(tmp,-param->powern) / (2.0*param->powern))/sqrt(tmp); + return (1.0 - pow(tmp,-1.0*param->powern) / (2.0*param->powern))/sqrt(tmp); if (tmp < param->c4) return 1.0; if (tmp < param->c3) return 1.0 - pow(tmp,param->powern)/(2.0*param->powern); @@ -1786,14 +1786,14 @@ double PairComb::qfo_self(Param *param, double qi, double selfpot) // sprintf(str,"Pair COMB charge %.10f with force %.10f hit min barrier", // qi,self_d); // error->warning(FLERR,str,0); - self_d += 4.0 * cmin * pow((qi-qmin),3); + self_d += 4.0 * cmin * pow((qi-qmin),3.0); } if (qi > qmax) { // char str[128]; // sprintf(str,"Pair COMB charge %.10f with force %.10f hit max barrier", // qi,self_d); // error->warning(FLERR,str,0); - self_d += 4.0 * cmax * pow((qi-qmax),3); + self_d += 4.0 * cmax * pow((qi-qmax),3.0); } return self_d; @@ -1896,9 +1896,9 @@ void PairComb::qfo_short(Param *param, int i, int j, double rsq, Asi = param->biga1 * exp(param->lam11*Di); Asj = param->biga2 * exp(param->lam12*Dj); Bsi = param->bigb1 * exp(param->lam21*Di)* - (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10))); + (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10.0))); Bsj = param->bigb2 * exp(param->lam22*Dj)* - (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10))); + (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10.0))); QUchi = (param->QU1-qi)*param->bD1; QUchj = (param->QU2-qj)*param->bD2; diff --git a/src/MANYBODY/pair_lcbop.cpp b/src/MANYBODY/pair_lcbop.cpp index a933e32037..882cf1686d 100644 --- a/src/MANYBODY/pair_lcbop.cpp +++ b/src/MANYBODY/pair_lcbop.cpp @@ -909,7 +909,7 @@ double PairLCBOP::gSpline( double x, double *dgdc ) { double PairLCBOP::hSpline( double x, double *dhdx ) { if( x < -d ) { double z = kappa*( x+d ); - double y = pow(z, 10); + double y = pow(z, 10.0); double w = pow( 1+y, -0.1 ); *dhdx = kappa*L*w/(1+y); return L*( 1 + z*w ); diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index bc450dbacc..6e1bdf6113 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -79,7 +79,7 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) : (2.0*MY_PI*gas_mass*force->mvv2e* force->boltz*reservoir_temperature)); sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e); - zz = exp(beta*chemical_potential)/(pow(lambda,3)); + zz = exp(beta*chemical_potential)/(pow(lambda,3.0)); // set defaults diff --git a/src/MOLECULE/angle_cosine_periodic.cpp b/src/MOLECULE/angle_cosine_periodic.cpp index 2743b38b1f..f3b4290613 100644 --- a/src/MOLECULE/angle_cosine_periodic.cpp +++ b/src/MOLECULE/angle_cosine_periodic.cpp @@ -138,8 +138,8 @@ void AngleCosinePeriodic::compute(int eflag, int vflag) un_2 = un_1; un_1 = un; } - tn = b_factor*pow(-1.0,m)*tn; - un = b_factor*pow(-1.0,m)*m*un; + tn = b_factor*pow(-1.0,(double)m)*tn; + un = b_factor*pow(-1.0,(double)m)*m*un; if (eflag) eangle = 2*k[type]*(1.0 - tn); @@ -286,5 +286,5 @@ double AngleCosinePeriodic::single(int type, int i1, int i2, int i3) if (c < -1.0) c = -1.0; c = cos(acos(c)*multiplicity[type]); - return k[type]*(1.0-b[type]*pow(-1.0,multiplicity[type])*c); + return k[type]*(1.0-b[type]*pow(-1.0,(double)multiplicity[type])*c); } diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp index 97a127b318..5858e48556 100644 --- a/src/MOLECULE/improper_umbrella.cpp +++ b/src/MOLECULE/improper_umbrella.cpp @@ -273,7 +273,7 @@ void ImproperUmbrella::coeff(int narg, char **arg) kw[i] = k_one; w0[i] = w_one/180.0 * MY_PI; if (w_one == 0) C[i] = 1.0; - else C[i] = kw[i]/(pow(sin(w0[i]),2)); + else C[i] = kw[i]/(pow(sin(w0[i]),2.0)); setflag[i] = 1; count++; } diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 3ad2c4e351..77302e6694 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -169,9 +169,9 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; force_kernel = r10inv*(pm->lj1*r2inv - pm->lj2)*r2inv * - pow(c,pm->ap); + pow(c,(double)pm->ap); force_angle = pm->ap * r10inv*(pm->lj3*r2inv - pm->lj4) * - pow(c,pm->ap-1)*s; + pow(c,(double)pm->ap-1.0)*s; eng_lj = r10inv*(pm->lj3*r2inv - pm->lj4); if (rsq > pm->cut_innersq) { @@ -185,7 +185,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) } if (eflag) { - evdwl = eng_lj * pow(c,pm->ap); + evdwl = eng_lj * pow(c,(double)pm->ap); evdwl *= factor_hb; ehbond += evdwl; } @@ -510,9 +510,9 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; force_kernel = r10inv*(pm->lj1*r2inv - pm->lj2)*r2inv * - pow(c,pm->ap); + pow(c,(double)pm->ap); force_angle = pm->ap * r10inv*(pm->lj3*r2inv - pm->lj4) * - pow(c,pm->ap-1)*s; + pow(c,(double)pm->ap-1.0)*s; // only lj part for now @@ -527,8 +527,8 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, eng_lj *= switch1; } - fforce += force_kernel*pow(c,pm->ap) + eng_lj*force_angle; - eng += eng_lj * pow(c,pm->ap) * factor_hb; + fforce += force_kernel*pow(c,(double)pm->ap) + eng_lj*force_angle; + eng += eng_lj * pow(c,(double)pm->ap) * factor_hb; } return eng; diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index e7f9c1e7c5..185dbc84e0 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -139,8 +139,8 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) r = sqrt(rsq); dr = r - pm->r0; dexp = exp(-pm->alpha * dr); - force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,pm->ap); - force_angle = pm->ap * eng_morse * pow(c,pm->ap-1)*s; + force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,(double)pm->ap); + force_angle = pm->ap * eng_morse * pow(c,pm->ap-1.0)*s; eng_morse = pm->d0 * (dexp*dexp - 2.0*dexp); if (rsq > pm->cut_innersq) { @@ -154,7 +154,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) } if (eflag) { - evdwl = eng_morse * pow(c,params[m].ap); + evdwl = eng_morse * pow(c,(double)params[m].ap); evdwl *= factor_hb; ehbond += evdwl; } @@ -415,8 +415,8 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, r = sqrt(rsq); dr = r - pm->r0; dexp = exp(-pm->alpha * dr); - force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,pm->ap); - force_angle = pm->ap * eng_morse * pow(c,pm->ap-1)*s; + force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,(double)pm->ap); + force_angle = pm->ap * eng_morse * pow(c,pm->ap-1.0)*s; eng_morse = pm->d0 * (dexp*dexp - 2.0*dexp); if (rsq > pm->cut_innersq) { @@ -429,8 +429,8 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, eng_morse *= switch1; } - eng += eng_morse * pow(c,params[m].ap)* factor_hb; - fforce += force_kernel*pow(c,pm->ap) + eng_morse*force_angle; + eng += eng_morse * pow(c,(double)params[m].ap)* factor_hb; + fforce += force_kernel*pow(c,(double)pm->ap) + eng_morse*force_angle; } return eng; diff --git a/src/REAX/pair_reax.cpp b/src/REAX/pair_reax.cpp index 3bdf165336..30f787d8d0 100644 --- a/src/REAX/pair_reax.cpp +++ b/src/REAX/pair_reax.cpp @@ -699,7 +699,7 @@ void PairREAX::taper_setup() double swb2,swa2,swb3,swa3,d1,d7; d1=swb-swa; - d7=pow(d1,7); + d7=pow(d1,7.0); swa2=swa*swa; swa3=swa2*swa; swb2=swb*swb; diff --git a/src/USER-EFF/fix_nvt_sllod_eff.cpp b/src/USER-EFF/fix_nvt_sllod_eff.cpp index adb3dd8221..17d58664db 100644 --- a/src/USER-EFF/fix_nvt_sllod_eff.cpp +++ b/src/USER-EFF/fix_nvt_sllod_eff.cpp @@ -124,7 +124,7 @@ void FixNVTSllodEff::nh_v_temp() temperature->restore_bias(i,v[i]); if (fabs(spin[i])==1) ervel[i] = ervel[i]*factor_eta - - dthalf*sqrt(pow(vdelu[0],2)+pow(vdelu[1],2)+pow(vdelu[2],2)); + dthalf*sqrt(vdelu[0]*vdelu[0]+vdelu[1]*vdelu[1]+vdelu[2]*vdelu[2]); } } } diff --git a/src/USER-EFF/pair_eff_inline.h b/src/USER-EFF/pair_eff_inline.h index 6f7705b193..6bd44fb7c2 100644 --- a/src/USER-EFF/pair_eff_inline.h +++ b/src/USER-EFF/pair_eff_inline.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -485,7 +485,7 @@ inline void PauliCoreElec(double rc, double re2, double *epauli, double *frc, (ssq + PAULI_CORE_C)) / (ssq + PAULI_CORE_C); dEdre2 = 2 * PAULI_CORE_A * PAULI_CORE_B * re2 * rcsq * exp(-PAULI_CORE_B * rcsq / - (ssq + PAULI_CORE_C)) / pow((PAULI_CORE_C + ssq),2); + (ssq + PAULI_CORE_C)) / ((PAULI_CORE_C + ssq)*(PAULI_CORE_C + ssq)); *epauli += E; *frc -= dEdrc; diff --git a/src/USER-EWALDN/ewald_n.cpp b/src/USER-EWALDN/ewald_n.cpp index f9f448d9e5..b6291cb5de 100644 --- a/src/USER-EWALDN/ewald_n.cpp +++ b/src/USER-EWALDN/ewald_n.cpp @@ -1205,7 +1205,7 @@ double EwaldN::NewtonSolve(double x, double Rc, bigint natoms, double vol, doubl double EwaldN::f(double x, double Rc, bigint natoms, double vol, double b2) { double a = Rc*x; - double f = (4.0*MY_PI*b2*pow(x,4)/vol/sqrt(natoms)*erfc(a) * + double f = (4.0*MY_PI*b2*pow(x,4.0)/vol/sqrt((double)natoms)*erfc(a) * (6.0*pow(a,-5.0) + 6.0*pow(a,-3.0) + 3.0/a + a) - accuracy); return f; } diff --git a/src/USER-MISC/fix_imd.cpp b/src/USER-MISC/fix_imd.cpp index 4f02a6f5df..48606876f0 100644 --- a/src/USER-MISC/fix_imd.cpp +++ b/src/USER-MISC/fix_imd.cpp @@ -47,6 +47,16 @@ negotiate an appropriate license for such distribution." Justin Gullingsrud, and James Phillips, (TCBG, Beckman Institute, UIUC) ------------------------------------------------------------------------- */ +#include "fix_imd.h" +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "error.h" +#include "group.h" +#include "memory.h" + #include #include #include @@ -68,16 +78,6 @@ negotiate an appropriate license for such distribution." #include -#include "fix_imd.h" -#include "atom.h" -#include "comm.h" -#include "update.h" -#include "respa.h" -#include "domain.h" -#include "error.h" -#include "group.h" -#include "memory.h" - /* re-usable integer hash table code with static linkage. */ /** hash table top level data structure */ diff --git a/src/USER-MISC/improper_ring.cpp b/src/USER-MISC/improper_ring.cpp index ef6940cc2d..cead4c87df 100644 --- a/src/USER-MISC/improper_ring.cpp +++ b/src/USER-MISC/improper_ring.cpp @@ -174,7 +174,7 @@ void ImproperRing::compute(int eflag, int vflag) /* Append the current angle to the sum of angle differences. */ angle_summer += (bend_angle[icomb] - chi[type]); } - if (eflag) eimproper = (1.0/6.0) *k[type] * pow(angle_summer,6); + if (eflag) eimproper = (1.0/6.0) *k[type] * pow(angle_summer,6.0); /* printf("The tags: %d-%d-%d-%d, of type %d .\n",atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4],type); // printf("The coordinates of the first: %f, %f, %f.\n", x[i1][0], x[i1][1], x[i1][2]); @@ -188,7 +188,7 @@ void ImproperRing::compute(int eflag, int vflag) /* Force calculation acting on all atoms. Calculate the derivatives of the potential. */ - angfac = k[type] * pow(angle_summer,5); + angfac = k[type] * pow(angle_summer,5.0); f1[0] = 0.0; f1[1] = 0.0; f1[2] = 0.0; f3[0] = 0.0; f3[1] = 0.0; f3[2] = 0.0; diff --git a/src/USER-MISC/pair_dipole_sf.cpp b/src/USER-MISC/pair_dipole_sf.cpp index 8adf066352..bd9e900c52 100644 --- a/src/USER-MISC/pair_dipole_sf.cpp +++ b/src/USER-MISC/pair_dipole_sf.cpp @@ -260,7 +260,7 @@ void PairDipoleSF::compute(int eflag, int vflag) if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { ecoul = qtmp * q[j] * rinv * - pow((1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])),2); + pow((1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])),2.0); if (mu[i][3] > 0.0 && mu[j][3] > 0.0) ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr); if (mu[i][3] > 0.0 && q[j] != 0.0) diff --git a/src/USER-OMP/angle_cosine_periodic_omp.cpp b/src/USER-OMP/angle_cosine_periodic_omp.cpp index b96f456f69..0f121efd73 100644 --- a/src/USER-OMP/angle_cosine_periodic_omp.cpp +++ b/src/USER-OMP/angle_cosine_periodic_omp.cpp @@ -163,8 +163,8 @@ void AngleCosinePeriodicOMP::eval(int nfrom, int nto, ThrData * const thr) un_2 = un_1; un_1 = un; } - tn = b_factor*pow(-1.0,m)*tn; - un = b_factor*pow(-1.0,m)*m*un; + tn = b_factor*pow(-1.0,(double)m)*tn; + un = b_factor*pow(-1.0,(double)m)*m*un; if (EFLAG) eangle = 2*k[type]*(1.0 - tn); diff --git a/src/USER-OMP/improper_ring_omp.cpp b/src/USER-OMP/improper_ring_omp.cpp index 7afc6aeba4..99ae6e532b 100644 --- a/src/USER-OMP/improper_ring_omp.cpp +++ b/src/USER-OMP/improper_ring_omp.cpp @@ -172,7 +172,7 @@ void ImproperRingOMP::eval(int nfrom, int nto, ThrData * const thr) /* Append the current angle to the sum of angle differences. */ angle_summer += (bend_angle[icomb] - chi[type]); } - if (EFLAG) eimproper = (1.0/6.0) *k[type] * pow(angle_summer,6); + if (EFLAG) eimproper = (1.0/6.0) *k[type] * pow(angle_summer,6.0); /* printf("The tags: %d-%d-%d-%d, of type %d .\n",atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4],type); // printf("The coordinates of the first: %f, %f, %f.\n", x[i1][0], x[i1][1], x[i1][2]); @@ -186,7 +186,7 @@ void ImproperRingOMP::eval(int nfrom, int nto, ThrData * const thr) /* Force calculation acting on all atoms. Calculate the derivatives of the potential. */ - angfac = k[type] * pow(angle_summer,5); + angfac = k[type] * pow(angle_summer,5.0); f1[0] = 0.0; f1[1] = 0.0; f1[2] = 0.0; f3[0] = 0.0; f3[1] = 0.0; f3[2] = 0.0; diff --git a/src/USER-OMP/pair_beck_omp.cpp b/src/USER-OMP/pair_beck_omp.cpp index fe279e6b96..231aad0390 100644 --- a/src/USER-OMP/pair_beck_omp.cpp +++ b/src/USER-OMP/pair_beck_omp.cpp @@ -126,7 +126,7 @@ void PairBeckOMP::eval(int iifrom, int iito, ThrData * const thr) alphaij = alpha[itype][jtype]; betaij = beta[itype][jtype]; term1 = aaij*aaij + rsq; - term2 = 1.0/pow(term1,5); + term2 = 1.0/pow(term1,5.0); term3 = 21.672 + 30.0*aaij*aaij + 6.0*rsq; term4 = alphaij + r5*betaij; term5 = alphaij + 6.0*r5*betaij; @@ -146,7 +146,7 @@ void PairBeckOMP::eval(int iifrom, int iito, ThrData * const thr) } if (EFLAG) { - term6 = 1.0/pow(term1,3); + term6 = 1.0/pow(term1,3.0); term1inv = 1.0/term1; evdwl = AA[itype][jtype]*exp(-1.0*r*term4); evdwl -= BB[itype][jtype]*term6*(1.0+(2.709+3.0*aaij*aaij)*term1inv); diff --git a/src/USER-OMP/pair_brownian_omp.cpp b/src/USER-OMP/pair_brownian_omp.cpp index 7a719e8b08..8da8221345 100644 --- a/src/USER-OMP/pair_brownian_omp.cpp +++ b/src/USER-OMP/pair_brownian_omp.cpp @@ -105,11 +105,11 @@ void PairBrownianOMP::compute(int eflag, int vflag) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); + RT0 = 8*MY_PI*mu*pow(rad,3.0); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -254,7 +254,7 @@ void PairBrownianOMP::eval(int iifrom, int iito, ThrData * const thr) if (FLAGLOG) { a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); diff --git a/src/USER-OMP/pair_brownian_poly_omp.cpp b/src/USER-OMP/pair_brownian_poly_omp.cpp index c0438f6497..891d868c25 100644 --- a/src/USER-OMP/pair_brownian_poly_omp.cpp +++ b/src/USER-OMP/pair_brownian_poly_omp.cpp @@ -105,11 +105,11 @@ void PairBrownianPolyOMP::compute(int eflag, int vflag) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); + RT0 = 8*MY_PI*mu*pow(rad,3.0); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -250,20 +250,20 @@ void PairBrownianPolyOMP::eval(int iifrom, int iito, ThrData * const thr) if (FLAGLOG) { a_sq = beta0*beta0/beta1/beta1/h_sep + - (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); - a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) + - pow(beta0,4))/21.0/pow(beta1,4)*h_sep*log(1.0/h_sep); + (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3.0)*log(1.0/h_sep); + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3.0) + + pow(beta0,4.0))/21.0/pow(beta1,4.0)*h_sep*log(1.0/h_sep); a_sq *= 6.0*MY_PI*mu*radi; - a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) * + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3.0) * log(1.0/h_sep); - a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + - 16.0*pow(beta0,4))/375.0/pow(beta1,4) * + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3.0) + + 16.0*pow(beta0,4.0))/375.0/pow(beta1,4.0) * h_sep*log(1.0/h_sep); a_sh *= 6.0*MY_PI*mu*radi; a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 * - pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); - a_pu *= 8.0*MY_PI*mu*pow(radi,3); + pow(beta0,3.0))/250.0/pow(beta1,3.0)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3.0); } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); diff --git a/src/USER-OMP/pair_dipole_sf_omp.cpp b/src/USER-OMP/pair_dipole_sf_omp.cpp index 8ea425fc8a..07721562b8 100644 --- a/src/USER-OMP/pair_dipole_sf_omp.cpp +++ b/src/USER-OMP/pair_dipole_sf_omp.cpp @@ -275,7 +275,7 @@ void PairDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) if (EFLAG) { if (rsq < cut_coulsq[itype][jtype]) { ecoul = qtmp * q[j] * rinv * - pow((1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])),2); + pow((1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])),2.0); if (mu[i][3] > 0.0 && mu[j][3] > 0.0) ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr); if (mu[i][3] > 0.0 && q[j] != 0.0) diff --git a/src/USER-OMP/pair_hbond_dreiding_lj_omp.cpp b/src/USER-OMP/pair_hbond_dreiding_lj_omp.cpp index fc4f7eba0e..074271a193 100644 --- a/src/USER-OMP/pair_hbond_dreiding_lj_omp.cpp +++ b/src/USER-OMP/pair_hbond_dreiding_lj_omp.cpp @@ -212,9 +212,9 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; force_kernel = r10inv*(pm->lj1*r2inv - pm->lj2)*r2inv * - pow(c,pm->ap); + pow(c,(double)pm->ap); force_angle = pm->ap * r10inv*(pm->lj3*r2inv - pm->lj4) * - pow(c,pm->ap-1)*s; + pow(c,pm->ap-1.0)*s; eng_lj = r10inv*(pm->lj3*r2inv - pm->lj4); if (rsq > pm->cut_innersq) { @@ -228,7 +228,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) } if (EFLAG) { - evdwl = eng_lj * pow(c,pm->ap); + evdwl = eng_lj * pow(c,(double)pm->ap); evdwl *= factor_hb; } diff --git a/src/USER-OMP/pair_hbond_dreiding_morse_omp.cpp b/src/USER-OMP/pair_hbond_dreiding_morse_omp.cpp index 291825ea78..fe0106473f 100644 --- a/src/USER-OMP/pair_hbond_dreiding_morse_omp.cpp +++ b/src/USER-OMP/pair_hbond_dreiding_morse_omp.cpp @@ -211,8 +211,8 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) r = sqrt(rsq); dr = r - pm->r0; dexp = exp(-pm->alpha * dr); - force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,pm->ap); - force_angle = pm->ap * eng_morse * pow(c,pm->ap-1)*s; + force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,(double)pm->ap); + force_angle = pm->ap * eng_morse * pow(c,(double)pm->ap-1.0)*s; eng_morse = pm->d0 * (dexp*dexp - 2.0*dexp); if (rsq > pm->cut_innersq) { @@ -226,7 +226,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) } if (EFLAG) { - evdwl = eng_morse * pow(c,params[m].ap); + evdwl = eng_morse * pow(c,(double)params[m].ap); evdwl *= factor_hb; } diff --git a/src/USER-OMP/pair_lubricate_omp.cpp b/src/USER-OMP/pair_lubricate_omp.cpp index 8808314c53..00a81f345b 100644 --- a/src/USER-OMP/pair_lubricate_omp.cpp +++ b/src/USER-OMP/pair_lubricate_omp.cpp @@ -99,13 +99,13 @@ void PairLubricateOMP::compute(int eflag, int vflag) double vol_f = vol_P/vol_T; if (flaglog == 0) { R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)* + RT0 = 8*MY_PI*mu*pow(rad,3.0); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)* (1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); - RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)* + RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)* (1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } @@ -344,7 +344,7 @@ void PairLubricateOMP::eval(int iifrom, int iito, ThrData * const thr) if (FLAGLOG) { a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3.0)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); diff --git a/src/USER-OMP/pair_lubricate_poly_omp.cpp b/src/USER-OMP/pair_lubricate_poly_omp.cpp index 20bb81af52..306c3f9bdf 100644 --- a/src/USER-OMP/pair_lubricate_poly_omp.cpp +++ b/src/USER-OMP/pair_lubricate_poly_omp.cpp @@ -351,21 +351,21 @@ void PairLubricatePolyOMP::eval(int iifrom, int iito, ThrData * const thr) if (FLAGLOG) { a_sq = beta0*beta0/beta1/beta1/h_sep + - (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); + (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3.0)*log(1.0/h_sep); a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0 * - pow(beta0,3)+pow(beta0,4))/21.0/pow(beta1,4) * + pow(beta0,3.0)+pow(beta0,4.0))/21.0/pow(beta1,4.0) * h_sep*log(1.0/h_sep); a_sq *= 6.0*MY_PI*mu*radi; - a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) * + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3.0) * log(1.0/h_sep); - a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + - 16.0*pow(beta0,4))/375.0/pow(beta1,4) * + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3.0) + + 16.0*pow(beta0,4.0))/375.0/pow(beta1,4.0) * h_sep*log(1.0/h_sep); a_sh *= 6.0*MY_PI*mu*radi; a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 * - pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); - a_pu *= 8.0*MY_PI*mu*pow(radi,3); + pow(beta0,3.0))/250.0/pow(beta1,3.0)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3.0); } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); // relative velocity at the point of closest approach diff --git a/src/USER-OMP/pppm_cg_omp.cpp b/src/USER-OMP/pppm_cg_omp.cpp index aa15a96517..47708f9829 100644 --- a/src/USER-OMP/pppm_cg_omp.cpp +++ b/src/USER-OMP/pppm_cg_omp.cpp @@ -264,14 +264,15 @@ void PPPMCGOMP::setup() if (sqk != 0.0) { numerator = form*MY_4PI/sqk; - denominator = gf_denom(snx,sny,snz); + denominator = gf_denom(snx,sny,snz); + const double dorder = static_cast(order); sum1 = 0.0; for (nx = -nbx; nx <= nbx; nx++) { qx = fkxk + unitkx*nx_pppm*nx; sx = exp(qx*qx/gew2); wx = 1.0; argx = 0.5*qx*xprd/nx_pppm; - if (argx != 0.0) wx = pow(sin(argx)/argx,order); + if (argx != 0.0) wx = pow(sin(argx)/argx,dorder); wx *=wx; for (ny = -nby; ny <= nby; ny++) { @@ -279,7 +280,7 @@ void PPPMCGOMP::setup() sy = exp(qy*qy/gew2); wy = 1.0; argy = 0.5*qy*yprd/ny_pppm; - if (argy != 0.0) wy = pow(sin(argy)/argy,order); + if (argy != 0.0) wy = pow(sin(argy)/argy,dorder); wy *= wy; for (nz = -nbz; nz <= nbz; nz++) { @@ -287,7 +288,7 @@ void PPPMCGOMP::setup() sz = exp(qz*qz/gew2); wz = 1.0; argz = 0.5*qz*zprd_slab/nz_pppm; - if (argz != 0.0) wz = pow(sin(argz)/argz,order); + if (argz != 0.0) wz = pow(sin(argz)/argz,dorder); wz *= wz; dot1 = fkxk*qx + fkyl*qy + fkzm*qz; diff --git a/src/USER-OMP/pppm_omp.cpp b/src/USER-OMP/pppm_omp.cpp index 496ee27aa8..24308de534 100644 --- a/src/USER-OMP/pppm_omp.cpp +++ b/src/USER-OMP/pppm_omp.cpp @@ -267,12 +267,13 @@ void PPPMOMP::setup() numerator = form*MY_4PI/sqk; denominator = gf_denom(snx,sny,snz); sum1 = 0.0; + const double dorder = static_cast(order); for (nx = -nbx; nx <= nbx; nx++) { qx = fkxk + unitkx*nx_pppm*nx; sx = exp(qx*qx/gew2); wx = 1.0; argx = 0.5*qx*xprd/nx_pppm; - if (argx != 0.0) wx = pow(sin(argx)/argx,order); + if (argx != 0.0) wx = pow(sin(argx)/argx,dorder); wx *=wx; for (ny = -nby; ny <= nby; ny++) { @@ -280,7 +281,7 @@ void PPPMOMP::setup() sy = exp(qy*qy/gew2); wy = 1.0; argy = 0.5*qy*yprd/ny_pppm; - if (argy != 0.0) wy = pow(sin(argy)/argy,order); + if (argy != 0.0) wy = pow(sin(argy)/argy,dorder); wy *= wy; for (nz = -nbz; nz <= nbz; nz++) { @@ -288,7 +289,7 @@ void PPPMOMP::setup() sz = exp(qz*qz/gew2); wz = 1.0; argz = 0.5*qz*zprd_slab/nz_pppm; - if (argz != 0.0) wz = pow(sin(argz)/argz,order); + if (argz != 0.0) wz = pow(sin(argz)/argz,dorder); wz *= wz; dot1 = fkxk*qx + fkyl*qy + fkzm*qz; diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index 21de2c0e1d..cde118fa33 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -338,12 +338,11 @@ void FixQEqReax::init_taper() else if (swb < 5 && comm->me == 0) error->warning(FLERR,"Fix qeq/reax has very low Taper radius cutoff"); - d7 = pow( swb - swa, 7 ); + d7 = pow( swb - swa, 7.0 ); swa2 = SQR( swa ); swa3 = CUBE( swa ); swb2 = SQR( swb ); swb3 = CUBE( swb ); - Tap[7] = 20.0 / d7; Tap[6] = -70.0 * (swa + swb) / d7; Tap[5] = 84.0 * (swa2 + 3.0*swa*swb + swb2) / d7; diff --git a/src/USER-REAXC/reaxc_valence_angles.cpp b/src/USER-REAXC/reaxc_valence_angles.cpp index 47d562ac82..6eee1548a1 100644 --- a/src/USER-REAXC/reaxc_valence_angles.cpp +++ b/src/USER-REAXC/reaxc_valence_angles.cpp @@ -60,7 +60,7 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk, real sqr_d_ji = SQR(d_ji); real sqr_d_jk = SQR(d_jk); real inv_dists = 1.0 / (d_ji * d_jk); - real inv_dists3 = pow( inv_dists, 3 ); + real inv_dists3 = pow( inv_dists, 3.0 ); real dot_dvecs = Dot( dvec_ji, dvec_jk, 3 ); real Cdot_inv3 = dot_dvecs * inv_dists3; diff --git a/src/lmpwindows.h b/src/lmpwindows.h index e865385750..8cecf8a7e3 100644 --- a/src/lmpwindows.h +++ b/src/lmpwindows.h @@ -24,10 +24,6 @@ inline double pow(int i, int j){ return pow((double)i,(double) j); } - -inline double pow(double i, int j){ - return pow(i,(double) j); -} #else inline double pow(int i, int j){ return pow((double)i,j); diff --git a/src/pair.cpp b/src/pair.cpp index 7a1b9d1153..4cdc7a8542 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1151,19 +1151,19 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits, error->warning(FLERR,"Table inner cutoff >= outer cutoff"); int nlowermin = 1; - while (!((pow(double(2),nlowermin) <= inner*inner) && - (pow(double(2),nlowermin+1) > inner*inner))) { - if (pow(double(2),nlowermin) <= inner*inner) nlowermin++; + while (!((pow(double(2),(double)nlowermin) <= inner*inner) && + (pow(double(2),(double)nlowermin+1.0) > inner*inner))) { + if (pow(double(2),(double)nlowermin) <= inner*inner) nlowermin++; else nlowermin--; } int nexpbits = 0; - double required_range = outer*outer / pow(double(2),nlowermin); + double required_range = outer*outer / pow(double(2),(double)nlowermin); double available_range = 2.0; while (available_range < required_range) { nexpbits++; - available_range = pow(double(2),pow(double(2),nexpbits)); + available_range = pow(double(2),pow(double(2),(double)nexpbits)); } int nmantbits = ntablebits - nexpbits; diff --git a/src/pair_beck.cpp b/src/pair_beck.cpp index 443de0cccf..d96cae416c 100644 --- a/src/pair_beck.cpp +++ b/src/pair_beck.cpp @@ -105,7 +105,7 @@ void PairBeck::compute(int eflag, int vflag) alphaij = alpha[itype][jtype]; betaij = beta[itype][jtype]; term1 = aaij*aaij + rsq; - term2 = 1.0/pow(term1,5); + term2 = 1.0/pow(term1,5.0); term3 = 21.672 + 30.0*aaij*aaij + 6.0*rsq; term4 = alphaij + r5*betaij; term5 = alphaij + 6.0*r5*betaij; @@ -125,7 +125,7 @@ void PairBeck::compute(int eflag, int vflag) } if (eflag) { - term6 = 1.0/pow(term1,3); + term6 = 1.0/pow(term1,3.0); term1inv = 1.0/term1; evdwl = AA[itype][jtype]*exp(-1.0*r*term4); evdwl -= BB[itype][jtype]*term6*(1.0+(2.709+3.0*aaij*aaij)*term1inv); @@ -344,7 +344,7 @@ double PairBeck::single(int i, int j, int itype, int jtype, alphaij = alpha[itype][jtype]; betaij = beta[itype][jtype]; term1 = aaij*aaij + rsq; - term2 = 1.0/pow(term1,5); + term2 = 1.0/pow(term1,5.0); term3 = 21.672 + 30.0*aaij*aaij + 6.0*rsq; term4 = alphaij + r5*betaij; term5 = alphaij + 6.0*r5*betaij; @@ -353,7 +353,7 @@ double PairBeck::single(int i, int j, int itype, int jtype, force_beck -= BB[itype][jtype]*r*term2*term3; fforce = factor_lj*force_beck*rinv; - term6 = 1.0/pow(term1,3); + term6 = 1.0/pow(term1,3.0); term1inv = 1.0/term1; phi_beck = AA[itype][jtype]*exp(-1.0*r*term4); phi_beck -= BB[itype][jtype]*term6*(1.0+(2.709+3.0*aaij*aaij)*term1inv);