jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2012-06-01 19:38:10 +00:00
parent 611e614170
commit 4990f2e53d
39 changed files with 186 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/FLD/pair_brownian.cpp
View File

@ -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;

29
src/FLD/pair_brownian_poly.cpp
View File

@ -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;

20
src/FLD/pair_lubricate.cpp
View File

@ -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);
}

30
src/FLD/pair_lubricateU.cpp
View File

@ -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);
}
}

11
src/FLD/pair_lubricateU_poly.cpp
View File

@ -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;

27
src/FLD/pair_lubricate_poly.cpp
View File

@ -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;

17
src/KSPACE/pppm.cpp
View File

@ -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<double>(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<double>(nx_pppm);
h_y = yprd/static_cast<double>(ny_pppm);
h_z = zprd_slab/static_cast<double>(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;
}

2
src/MANYBODY/pair_airebo.cpp
View File

@ -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 {

26
src/MANYBODY/pair_comb.cpp
View File

@ -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;

2
src/MANYBODY/pair_lcbop.cpp
View File

@ -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 );

2
src/MC/fix_gcmc.cpp
View File

@ -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

6
src/MOLECULE/angle_cosine_periodic.cpp
View File

@ -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);
}

2
src/MOLECULE/improper_umbrella.cpp
View File

@ -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++;
}

14
src/MOLECULE/pair_hbond_dreiding_lj.cpp
View File

@ -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;

14
src/MOLECULE/pair_hbond_dreiding_morse.cpp
View File

@ -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;

2
src/REAX/pair_reax.cpp
View File

@ -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;

2
src/USER-EFF/fix_nvt_sllod_eff.cpp
View File

@ -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]);
}
}
}

4
src/USER-EFF/pair_eff_inline.h
View File

@ -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;

2
src/USER-EWALDN/ewald_n.cpp
View File

@ -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;
}

20
src/USER-MISC/fix_imd.cpp
View File

@ -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 <math.h>
#include <stdio.h>
#include <stdlib.h>
@ -68,16 +78,6 @@ negotiate an appropriate license for such distribution."
#include <errno.h>
#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 */

4
src/USER-MISC/improper_ring.cpp
View File

@ -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;

2
src/USER-MISC/pair_dipole_sf.cpp
View File

@ -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)

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

@ -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);

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

@ -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;

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

@ -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);

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

@ -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);

20
src/USER-OMP/pair_brownian_poly_omp.cpp
View File

@ -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);

2
src/USER-OMP/pair_dipole_sf_omp.cpp
View File

@ -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)

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

@ -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;
}

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

@ -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;
}

10
src/USER-OMP/pair_lubricate_omp.cpp
View File

@ -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);

14
src/USER-OMP/pair_lubricate_poly_omp.cpp
View File

@ -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

9
src/USER-OMP/pppm_cg_omp.cpp
View File

@ -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<double>(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;

7
src/USER-OMP/pppm_omp.cpp
View File

@ -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<double>(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;

3
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -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;

2
src/USER-REAXC/reaxc_valence_angles.cpp
View File

@ -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;

4
src/lmpwindows.h
View File

@ -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);

10
src/pair.cpp
View File

@ -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;

8
src/pair_beck.cpp
View File

@ -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);