ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2013-02-08 16:57:30 +00:00
parent 6676724eee
commit ebe08bf4ca
9 changed files with 155 additions and 125 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/COLLOID/pair_colloid.cpp
View File

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

14
src/FLD/pair_brownian.cpp
View File

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

22
src/FLD/pair_brownian_poly.cpp
View File

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

36
src/MANYBODY/pair_airebo.cpp
View File

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

8
src/MOLECULE/angle_cosine_periodic.cpp
View File

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

62
src/MOLECULE/pair_hbond_dreiding_lj.cpp
View File

@ -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 = &params[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 = &params[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;

64
src/MOLECULE/pair_hbond_dreiding_morse.cpp
View File

@ -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 = &params[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 = &params[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;

46
src/math_special.h
View File

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

10
src/pair_beck.cpp
View File

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