From ebe08bf4ca16138334999a68a88b23700c28a1c8 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 8 Feb 2013 16:57:30 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9405 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/COLLOID/pair_colloid.cpp | 18 +++--- src/FLD/pair_brownian.cpp | 14 +++-- src/FLD/pair_brownian_poly.cpp | 22 ++++---- src/MANYBODY/pair_airebo.cpp | 36 ++++++------ src/MOLECULE/angle_cosine_periodic.cpp | 8 ++- src/MOLECULE/pair_hbond_dreiding_lj.cpp | 62 ++++++++++----------- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 64 +++++++++++----------- src/math_special.h | 46 +++++++++++----- src/pair_beck.cpp | 10 ++-- 9 files changed, 155 insertions(+), 125 deletions(-) diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index fb0d44df42..5ec424d581 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -27,8 +27,10 @@ #include "neigh_list.h" #include "memory.h" #include "error.h" +#include "math_special.h" using namespace LAMMPS_NS; +using namespace MathSpecial; /* ---------------------------------------------------------------------- */ @@ -155,10 +157,10 @@ void PairColloid::compute(int eflag, int vflag) K[6] = K[2]-r; K[7] = 1.0/(K[3]*K[4]); K[8] = 1.0/(K[5]*K[6]); - g[0] = pow(K[3],-7.0); - g[1] = pow(K[4],-7.0); - g[2] = pow(K[5],-7.0); - g[3] = pow(K[6],-7.0); + g[0] = powint(K[3],-7); + g[1] = powint(K[4],-7); + g[2] = powint(K[5],-7); + g[3] = powint(K[6],-7); h[0] = ((K[3]+5.0*K[1])*K[3]+30.0*K[0])*g[0]; h[1] = ((K[4]+5.0*K[1])*K[4]+30.0*K[0])*g[1]; h[2] = ((K[5]+5.0*K[2])*K[5]-30.0*K[0])*g[2]; @@ -488,10 +490,10 @@ double PairColloid::single(int i, int j, int itype, int jtype, double rsq, K[6] = K[2]-r; K[7] = 1.0/(K[3]*K[4]); K[8] = 1.0/(K[5]*K[6]); - g[0] = pow(K[3],-7.0); - g[1] = pow(K[4],-7.0); - g[2] = pow(K[5],-7.0); - g[3] = pow(K[6],-7.0); + g[0] = powint(K[3],-7); + g[1] = powint(K[4],-7); + g[2] = powint(K[5],-7); + g[3] = powint(K[6],-7); h[0] = ((K[3]+5.0*K[1])*K[3]+30.0*K[0])*g[0]; h[1] = ((K[4]+5.0*K[1])*K[4]+30.0*K[0])*g[1]; h[2] = ((K[5]+5.0*K[2])*K[5]-30.0*K[0])*g[2]; diff --git a/src/FLD/pair_brownian.cpp b/src/FLD/pair_brownian.cpp index 350bb997d1..4f81110c20 100755 --- a/src/FLD/pair_brownian.cpp +++ b/src/FLD/pair_brownian.cpp @@ -37,11 +37,13 @@ #include "variable.h" #include "random_mars.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; // same as fix_wall.cpp @@ -129,11 +131,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.0); + RT0 = 8*MY_PI*mu*cube(rad); //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.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*cube(rad)*(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 +210,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.0)*(3.0/160.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*cube(radi)*(3.0/160.0*log(1.0/h_sep)); } else a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); @@ -543,7 +545,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.0); + vol_P = atom->natoms*(4.0/3.0)*MY_PI*cube(rad); double vol_f = vol_P/vol_T; @@ -552,10 +554,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.0); // not actually needed + RT0 = 8*MY_PI*mu*cube(rad); // 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.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); } } diff --git a/src/FLD/pair_brownian_poly.cpp b/src/FLD/pair_brownian_poly.cpp index 7373f15de1..43a0d80628 100644 --- a/src/FLD/pair_brownian_poly.cpp +++ b/src/FLD/pair_brownian_poly.cpp @@ -38,11 +38,13 @@ #include "variable.h" #include "random_mars.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; // same as fix_wall.cpp @@ -115,11 +117,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.0); + RT0 = 8*MY_PI*mu*cube(rad); //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.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RT0 = 8*MY_PI*mu*cube(rad)*(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); } } @@ -200,20 +202,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.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); + (1.0+7.0*beta0+beta0*beta0)/5.0/cube(beta1)*log(1.0/h_sep); + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*cube(beta0) + + powint(beta0,4))/21.0/powint(beta1,4)*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.0) * + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/cube(beta1) * log(1.0/h_sep); - 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) * + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*cube(beta0) + + 16.0*powint(beta0,4))/375.0/powint(beta1,4) * 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.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); + cube(beta0))/250.0/cube(beta1)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*cube(radi); } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index a39d2484c4..6e1edb1ab5 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -32,11 +32,13 @@ #include "neigh_list.h" #include "neigh_request.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; #define MAXLINE 1024 #define TOL 1.0e-9 @@ -994,7 +996,7 @@ void PairAIREBO::TORSION(int eflag, int vflag) cw2 = (.5*(1.0-cw)); ekijl = epsilonT[ktype][ltype]; Ec = 256.0*ekijl/405.0; - Vtors = (Ec*(pow(cw2,5.0)))-(ekijl/10.0); + Vtors = (Ec*(powint(cw2,5)))-(ekijl/10.0); if (eflag) evdwl = Vtors*w21*w23*w34*(1.0-tspjik)*(1.0-tspijl); @@ -1048,7 +1050,7 @@ void PairAIREBO::TORSION(int eflag, int vflag) ddndil = cross321mag*dxjdil; dcwddn = -cwnum/(cwnom*cwnom); dcwdn = 1.0/cwnom; - dvpdcw = (-1.0)*Ec*(-.5)*5.0*pow(cw2,4.0) * + dvpdcw = (-1.0)*Ec*(-.5)*5.0*powint(cw2,4) * w23*w21*w34*(1.0-tspjik)*(1.0-tspijl); Ftmp[0] = dvpdcw*((dcwdn*dndij[0])+(dcwddn*ddndij*del23[0]/r23)); @@ -1289,7 +1291,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], dN2[1] = 0.0; PijS = PijSpline(NijC,NijH,itype,jtype,dN2); pij = pow(1.0+Etmp+PijS,-0.5); - tmp = -0.5*pow(pij,3.0); + tmp = -0.5*cube(pij); // pij forces @@ -1434,7 +1436,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], dN2[1] = 0.0; PjiS = PijSpline(NjiC,NjiH,jtype,itype,dN2); pji = pow(1.0+Etmp+PjiS,-0.5); - tmp = -0.5*pow(pji,3.0); + tmp = -0.5*cube(pji); REBO_neighs = REBO_firstneigh[j]; for (l = 0; l < REBO_numneigh[j]; l++) { @@ -1460,11 +1462,11 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) - (cosijl*rij[2]/(rijmag*rijmag)); dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) + - (cosijl*((rij[0]/pow(rijmag,2.0))-(rjl[0]/(rjlmag*rjlmag)))); + (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag)))); dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) + - (cosijl*((rij[1]/pow(rijmag,2.0))-(rjl[1]/(rjlmag*rjlmag)))); + (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag)))); dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) + - (cosijl*((rij[2]/pow(rijmag,2.0))-(rjl[2]/(rjlmag*rjlmag)))); + (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag)))); dcosijldrl[0] = (rij[0]/(rijmag*rjlmag))+(cosijl*rjl[0]/(rjlmag*rjlmag)); dcosijldrl[1] = (rij[1]/(rijmag*rjlmag))+(cosijl*rjl[1]/(rjlmag*rjlmag)); dcosijldrl[2] = (rij[2]/(rijmag*rjlmag))+(cosijl*rjl[2]/(rjlmag*rjlmag)); @@ -1774,7 +1776,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], cwnom = r21mag*r34mag*r23mag*r23mag*sin321*sin234; om1234 = cwnum/cwnom; cw = om1234; - Etmp += ((1.0-pow(om1234,2.0))*w21*w34) * + Etmp += ((1.0-square(om1234))*w21*w34) * (1.0-tspjik)*(1.0-tspijl); dt1dik = (rik2i)-(dctik*sink2i*cos321); @@ -1801,7 +1803,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], aa = (prefactor*2.0*cw/cwnom)*w21*w34 * (1.0-tspjik)*(1.0-tspijl); - aaa1 = -prefactor*(1.0-pow(om1234,2.0)) * + aaa1 = -prefactor*(1.0-square(om1234)) * (1.0-tspjik)*(1.0-tspijl); aaa2 = aaa1*w21*w34; at2 = aa*cwnum; @@ -2116,7 +2118,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, dN2PIJ[1] = 0.0; PijS = PijSpline(NijC,NijH,itype,jtype,dN2PIJ); pij = pow(1.0+Etmp+PijS,-0.5); - tmppij = -.5*pow(pij,3.0); + tmppij = -.5*cube(pij); tmp3pij = tmp3; tmp = 0.0; tmp2 = 0.0; @@ -2156,7 +2158,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, dN2PJI[1] = 0.0; PjiS = PijSpline(NjiC,NjiH,jtype,itype,dN2PJI); pji = pow(1.0+Etmp+PjiS,-0.5); - tmppji = -.5*pow(pji,3.0); + tmppji = -.5*cube(pji); tmp3pji = tmp3; // evaluate Nij conj @@ -2238,7 +2240,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, (crosskij[1]*crossijl[1]) + (crosskij[2]*crossijl[2])) / (crosskijmag*crossijlmag)); - Etmp += ((1.0-pow(omkijl,2.0))*wik*wjl) * + Etmp += ((1.0-square(omkijl))*wik*wjl) * (1.0-tspjik)*(1.0-tspijl); } } @@ -2392,11 +2394,11 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) - (cosijl*rij[2]/(rijmag*rijmag)); dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) + - (cosijl*((rij[0]/pow(rijmag,2.0))-(rjl[0]/(rjlmag*rjlmag)))); + (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag)))); dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) + - (cosijl*((rij[1]/pow(rijmag,2.0))-(rjl[1]/(rjlmag*rjlmag)))); + (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag)))); dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) + - (cosijl*((rij[2]/pow(rijmag,2.0))-(rjl[2]/(rjlmag*rjlmag)))); + (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag)))); dcosijldrl[0] = (rij[0]/(rijmag*rjlmag)) + (cosijl*rjl[0]/(rjlmag*rjlmag)); dcosijldrl[1] = (rij[1]/(rijmag*rjlmag)) + @@ -2703,7 +2705,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, cwnom = r21mag*r34mag*r23mag*r23mag*sin321*sin234; om1234 = cwnum/cwnom; cw = om1234; - Etmp += ((1.0-pow(om1234,2.0))*w21*w34) * + Etmp += ((1.0-square(om1234))*w21*w34) * (1.0-tspjik)*(1.0-tspijl); dt1dik = (rik2i)-(dctik*sink2i*cos321); @@ -2733,7 +2735,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, aa = (prefactor*2.0*cw/cwnom)*w21*w34 * (1.0-tspjik)*(1.0-tspijl); - aaa1 = -prefactor*(1.0-pow(om1234,2.0)) * + aaa1 = -prefactor*(1.0-square(om1234)) * (1.0-tspjik)*(1.0-tspijl); aaa2 = aaa1*w21*w34; at2 = aa*cwnum; diff --git a/src/MOLECULE/angle_cosine_periodic.cpp b/src/MOLECULE/angle_cosine_periodic.cpp index 39ab78114f..2d10052e33 100644 --- a/src/MOLECULE/angle_cosine_periodic.cpp +++ b/src/MOLECULE/angle_cosine_periodic.cpp @@ -24,11 +24,13 @@ #include "comm.h" #include "force.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; #define SMALL 0.001 @@ -136,8 +138,8 @@ void AngleCosinePeriodic::compute(int eflag, int vflag) un_2 = un_1; un_1 = un; } - tn = b_factor*pow(-1.0,(double)m)*tn; - un = b_factor*pow(-1.0,(double)m)*m*un; + tn = b_factor*powsign(m)*tn; + un = b_factor*powsign(m)*m*un; if (eflag) eangle = 2*k[type]*(1.0 - tn); @@ -284,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,(double)multiplicity[type])*c); + return k[type]*(1.0-b[type]*powsign(multiplicity[type])*c); } diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 200b344e1d..4239bea8a4 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -28,11 +28,13 @@ #include "neigh_list.h" #include "domain.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; #define SMALL 0.001 #define CHUNK 8 @@ -83,7 +85,6 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) double r2inv,r10inv; double switch1,switch2; int *ilist,*jlist,*klist,*numneigh,**firstneigh; - Param *pm; evdwl = ehbond = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -135,9 +136,9 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) ktype = type[k]; m = type2param[itype][jtype][ktype]; if (m < 0) continue; - pm = ¶ms[m]; + const Param &pm = params[m]; - if (rsq < pm->cut_outersq) { + if (rsq < pm.cut_outersq) { delr1[0] = x[i][0] - x[k][0]; delr1[1] = x[i][1] - x[k][1]; delr1[2] = x[i][2] - x[k][2]; @@ -160,7 +161,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); - if (ac > pm->cut_angle && ac < (2.0*MY_PI - pm->cut_angle)) { + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -168,24 +169,24 @@ 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,(double)pm->ap); - force_angle = pm->ap * r10inv*(pm->lj3*r2inv - pm->lj4) * - pow(c,(double)pm->ap-1.0)*s; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * + powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * + powint(c,pm.ap-1)*s; - eng_lj = r10inv*(pm->lj3*r2inv - pm->lj4); - if (rsq > pm->cut_innersq) { - switch1 = (pm->cut_outersq-rsq) * (pm->cut_outersq-rsq) * - (pm->cut_outersq + 2.0*rsq - 3.0*pm->cut_innersq) / - pm->denom_vdw; - switch2 = 12.0*rsq * (pm->cut_outersq-rsq) * - (rsq-pm->cut_innersq) / pm->denom_vdw; + eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; force_kernel = force_kernel*switch1 + eng_lj*switch2; eng_lj *= switch1; } if (eflag) { - evdwl = eng_lj * pow(c,(double)pm->ap); + evdwl = eng_lj * powint(c,pm.ap); evdwl *= factor_hb; ehbond += evdwl; } @@ -450,7 +451,6 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, double switch1,switch2; double delr1[3],delr2[3]; int *klist; - Param *pm; double **x = atom->x; int **special = atom->special; @@ -477,7 +477,7 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, ktype = type[k]; m = type2param[itype][jtype][ktype]; if (m < 0) continue; - pm = ¶ms[m]; + const Param &pm = params[m]; delr1[0] = x[i][0] - x[k][0]; delr1[1] = x[i][1] - x[k][1]; @@ -501,7 +501,7 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); - if (ac < pm->cut_angle || ac > (2.0*MY_PI - pm->cut_angle)) return 0.0; + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -509,26 +509,24 @@ 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,(double)pm->ap); - force_angle = pm->ap * r10inv*(pm->lj3*r2inv - pm->lj4) * - pow(c,(double)pm->ap-1.0)*s; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * + powint(c,pm.ap-1)*s; // only lj part for now - eng_lj = r10inv*(pm->lj3*r2inv - pm->lj4); - if (rsq > pm->cut_innersq) { - switch1 = (pm->cut_outersq-rsq) * (pm->cut_outersq-rsq) * - (pm->cut_outersq + 2.0*rsq - 3.0*pm->cut_innersq) / - pm->denom_vdw; - switch2 = 12.0*rsq * (pm->cut_outersq-rsq) * - (rsq-pm->cut_innersq) / pm->denom_vdw; + eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; force_kernel = force_kernel*switch1 + eng_lj*switch2; eng_lj *= switch1; } - fforce += force_kernel*pow(c,(double)pm->ap) + eng_lj*force_angle; - eng += eng_lj * pow(c,(double)pm->ap) * factor_hb; + fforce += force_kernel*powint(c,pm.ap) + eng_lj*force_angle; + eng += eng_lj * powint(c,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 07010c588d..d7f0f25993 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -28,11 +28,13 @@ #include "neigh_list.h" #include "domain.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; #define SMALL 0.001 #define CHUNK 8 @@ -53,7 +55,6 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) double fi[3],fj[3],delr1[3],delr2[3]; double r,dr,dexp,eng_morse,switch1,switch2; int *ilist,*jlist,*klist,*numneigh,**firstneigh; - Param *pm; evdwl = ehbond = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -105,9 +106,9 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) ktype = type[k]; m = type2param[itype][jtype][ktype]; if (m < 0) continue; - pm = ¶ms[m]; + const Param &pm = params[m]; - if (rsq < pm->cut_outersq) { + if (rsq < pm.cut_outersq) { delr1[0] = x[i][0] - x[k][0]; delr1[1] = x[i][1] - x[k][1]; delr1[2] = x[i][2] - x[k][2]; @@ -130,31 +131,31 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); - if (ac > pm->cut_angle && ac < (2.0*MY_PI - pm->cut_angle)) { + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; // Morse-specific kernel r = sqrt(rsq); - dr = r - pm->r0; - dexp = exp(-pm->alpha * dr); - eng_morse = pm->d0 * (dexp*dexp - 2.0*dexp); - 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; + dr = r - pm.r0; + dexp = exp(-pm.alpha * dr); + eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); + force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; - if (rsq > pm->cut_innersq) { - switch1 = (pm->cut_outersq-rsq) * (pm->cut_outersq-rsq) * - (pm->cut_outersq + 2.0*rsq - 3.0*pm->cut_innersq) / - pm->denom_vdw; - switch2 = 12.0*rsq * (pm->cut_outersq-rsq) * - (rsq-pm->cut_innersq) / pm->denom_vdw; + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; force_kernel = force_kernel*switch1 + eng_morse*switch2; eng_morse *= switch1; } if (eflag) { - evdwl = eng_morse * pow(c,(double)params[m].ap); + evdwl = eng_morse * powint(c,pm.ap); evdwl *= factor_hb; ehbond += evdwl; } @@ -355,7 +356,6 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, double switch1,switch2; double delr1[3],delr2[3]; int *klist; - Param *pm; double **x = atom->x; int **special = atom->special; @@ -382,7 +382,7 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, ktype = type[k]; m = type2param[itype][jtype][ktype]; if (m < 0) continue; - pm = ¶ms[m]; + const Param &pm = params[m]; delr1[0] = x[i][0] - x[k][0]; delr1[1] = x[i][1] - x[k][1]; @@ -406,31 +406,31 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); - if (ac < pm->cut_angle || ac > (2.0*MY_PI - pm->cut_angle)) return 0.0; + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; // Morse-specific kernel r = sqrt(rsq); - dr = r - pm->r0; - dexp = exp(-pm->alpha * dr); - eng_morse = pm->d0 * (dexp*dexp - 2.0*dexp); //<-- BUGFIX 2012-11-14 - 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; + dr = r - pm.r0; + dexp = exp(-pm.alpha * dr); + eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); //<-- BUGFIX 2012-11-14 + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); + force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; - if (rsq > pm->cut_innersq) { - switch1 = (pm->cut_outersq-rsq) * (pm->cut_outersq-rsq) * - (pm->cut_outersq + 2.0*rsq - 3.0*pm->cut_innersq) / - pm->denom_vdw; - switch2 = 12.0*rsq * (pm->cut_outersq-rsq) * - (rsq-pm->cut_innersq) / pm->denom_vdw; + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; force_kernel = force_kernel*switch1 + eng_morse*switch2; eng_morse *= switch1; } - eng += eng_morse * pow(c,(double)params[m].ap)* factor_hb; - fforce += force_kernel*pow(c,(double)pm->ap) + eng_morse*force_angle; + eng += eng_morse * powint(c,pm.ap)* factor_hb; + fforce += force_kernel*powint(c,pm.ap) + eng_morse*force_angle; } return eng; diff --git a/src/math_special.h b/src/math_special.h index 1062d9e7a5..c2fee94f2b 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -19,25 +19,45 @@ namespace LAMMPS_NS { namespace MathSpecial { - // x**2; + + // x**2, use instead of pow(x,2.0) + static inline double square(const double &x) { return x*x; } - // optimized version of (sin(x)/x)**n with n being a positive integer - static inline double powsinxx(const double x, int n) { - double xx,yy,ww; + // x**3, use instead of pow(x,3.0) + static inline double cube(const double &x) { return x*x*x; } - if ((x == 0.0) || (n == 0)) return 1.0; + // return -1.0 for odd n, 1.0 for even n, like pow(-1.0,n) + static inline double powsign(const int n) { return (n & 1) ? -1.0 : 1.0; } - xx = sin(x)/x; - yy = (n & 1) ? xx : 1.0; - ww = xx; - n >>= 1; + // optimized version of pow(x,n) with n being integer + // up to 10x faster than pow(x,y) - while (n) { - ww *= ww; + static inline double powint(const double &x, const int n) { + double yy,ww; + + if (x == 0.0) return 0.0; + int nn = (n > 0) ? n : -n; + ww = x; + + for (yy = 1.0; nn != 0; nn >>= 1, ww *=ww) + if (nn & 1) yy *= ww; + + return (n > 0) ? yy : 1.0/yy; + } + + // optimized version of (sin(x)/x)**n with n being a _positive_ integer + + static inline double powsinxx(const double &x, int n) { + double yy,ww; + + if (x == 0.0) return 1.0; + + ww = sin(x)/x; + + for (yy = 1.0; n != 0; n >>= 1, ww *=ww) if (n & 1) yy *= ww; - n >>=1; - } + return yy; } } diff --git a/src/pair_beck.cpp b/src/pair_beck.cpp index b5aee0f6c2..d41b7e6931 100644 --- a/src/pair_beck.cpp +++ b/src/pair_beck.cpp @@ -25,8 +25,10 @@ #include "neigh_list.h" #include "memory.h" #include "error.h" +#include "math_special.h" using namespace LAMMPS_NS; +using namespace MathSpecial; /* ---------------------------------------------------------------------- */ @@ -105,7 +107,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.0); + term2 = powint(term1,-5); term3 = 21.672 + 30.0*aaij*aaij + 6.0*rsq; term4 = alphaij + r5*betaij; term5 = alphaij + 6.0*r5*betaij; @@ -125,7 +127,7 @@ void PairBeck::compute(int eflag, int vflag) } if (eflag) { - term6 = 1.0/pow(term1,3.0); + term6 = powint(term1,-3); 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); @@ -343,7 +345,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.0); + term2 = powint(term1,-5); term3 = 21.672 + 30.0*aaij*aaij + 6.0*rsq; term4 = alphaij + r5*betaij; term5 = alphaij + 6.0*r5*betaij; @@ -352,7 +354,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.0); + term6 = powint(term1,-3); 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);