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

about 1.5x speedup for pair style comb3 by using MathSpecial::powint()

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-09-16 00:13:28 -04:00
parent 2b6ff442d8
commit 7aa6241db5
1 changed files with 34 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
src/MANYBODY/pair_comb3.cpp
View File

@ -27,6 +27,7 @@
#include "group.h"
#include "math_const.h"
#include "math_extra.h"
#include "math_special.h"
#include "memory.h"
#include "my_page.h"
#include "neigh_list.h"
@ -42,6 +43,7 @@
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathExtra;
using namespace MathSpecial;
#define DELTA 4
#define PGDELTA 1
@ -666,7 +668,7 @@ void PairComb3::setup_params()
params[m].Qo = (params[m].QU+params[m].QL)/2.0; // (A22)
params[m].dQ = (params[m].QU-params[m].QL)/2.0; // (A21)
params[m].aB = 1.0 /
(1.0-pow(fabs(params[m].Qo/params[m].dQ),10)); // (A20)
(1.0-powint(fabs(params[m].Qo/params[m].dQ),10)); // (A20)
params[m].bB = pow(fabs(params[m].aB),0.1)/params[m].dQ; // (A19)
params[m].nD = log(params[m].DU/(params[m].DU-params[m].DL))/
log(params[m].QU/(params[m].QU-params[m].QL));
@ -1432,7 +1434,7 @@ void PairComb3::repulsive(Param *parami, Param *paramj, double rsq,
vrcs = 1.0; fvrcs = 0.0;
if (romi != 0.0 && r < addr) {
vrcs += romi * pow((1.0-r/addr),2.0);
vrcs += romi * square(1.0-r/addr);
fvrcs = romi * 2.0 * (r/addr-1.0)/addr;
fforce = fforce*vrcs - caj * tmp_fc * vrcs * fvrcs;
}
@ -1455,7 +1457,7 @@ double PairComb3::zeta(Param *parami, Param *paramj, double rsqij,
costheta = dot3(delrij,delrik) / (rij*rik);
rlm3 = parami->beta;
arg = pow(rlm3*(rij-rik),int(parami->powermint));
arg = powint(rlm3*(rij-rik),int(parami->powermint));
if (arg > 69.0776) ex_delr = 1.e30;
else if (arg < -69.0776) ex_delr = 0.0;
else ex_delr = exp(arg);
@ -1776,8 +1778,8 @@ double PairComb3::self(Param *param, double qi)
self_tmp = qi*(s1+qi*(s2+qi*(s3+qi*s4)));
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 * powint((qi-qmin),4);
if (qi > qmax) self_tmp += cmax * powint((qi-qmax),4);
return self_tmp;
}
@ -1802,8 +1804,8 @@ void PairComb3::comb_fa(double r, Param *parami, Param *paramj, double iq,
Di = Dj = Bsi = 0.0;
Di = parami->DU + pow(fabs(parami->bD*(parami->QU-qi)),parami->nD);
Dj = paramj->DU + pow(fabs(paramj->bD*(paramj->QU-qj)),paramj->nD);
YYBn = (parami->aB-fabs(pow(parami->bB*(qi-parami->Qo),10)));
YYBj = (paramj->aB-fabs(pow(paramj->bB*(qj-paramj->Qo),10)));
YYBn = (parami->aB-fabs(powint(parami->bB*(qi-parami->Qo),10)));
YYBj = (paramj->aB-fabs(powint(paramj->bB*(qj-paramj->Qo),10)));
if (YYBn*YYBj > 0.0) {
AlfDiAlfDj = exp(0.5*(parami->alfi*Di+paramj->alfi*Dj));
@ -2124,7 +2126,7 @@ void PairComb3::comb_zetaterm_d(double prefac_ij1, double prefac_ij2,
fc_k = comb_fc(rik,paramk);
dfc_k = comb_fc_d(rik,paramk);
rlm3 = parami->beta;
tmp = pow(rlm3*(rij-rik),mint);
tmp = powint(rlm3*(rij-rik),mint);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
@ -2155,7 +2157,7 @@ void PairComb3::comb_zetaterm_d(double prefac_ij1, double prefac_ij2,
com3k = 0.0;
}
ex_delr_d = mint*pow(rlm3,mint)*pow((rij-rik),(mint-1))*ex_delr; // com3
ex_delr_d = mint*powint(rlm3,mint)*powint((rij-rik),(mint-1))*ex_delr; // com3
scale3(-dfc_k*gijk*ex_delr,rik_hat,dri); // com1
scaleadd3(fc_k*gijk_d*ex_delr,dcosdri,dri,dri); // com2
scaleadd3(fc_k*gijk*ex_delr_d,rik_hat,dri,dri); // com3 cont'd
@ -2690,10 +2692,10 @@ void PairComb3::field(Param *parami, Param *paramj, double rsq, double iq,
pcmi1 = parami->pcmn1;
pcmi2 = parami->pcmn2;
rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));
rf3i = r3/(square(r3)+cube(pcmi1));
rcf3i = rc3/(square(rc3)+cube(pcmi1));
rf5i = r5/(square(r5)+powint(pcmi2,5));
rcf5i = rc5/(square(rc5)+powint(pcmi2,5));
drf3i = 3/r*rf3i-6*rsq*rf3i*rf3i;
drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
@ -2742,7 +2744,7 @@ void PairComb3::rad_calc(double r, Param *parami, Param *paramj,
vrad = pradx = prady = pradz = 0.0;
xrad = -comb_fc(r,parami)*parami->pcross + xcn;
yrad = -comb_fc(r,paramj)*paramj->pcross + ycn;
zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
zcon = 1.0 + square(kconjug) + square(lconjug);
if (xrad < 0.0) xrad = 0.0;
if (yrad < 0.0) yrad = 0.0;
@ -2893,14 +2895,14 @@ double PairComb3::bbtor1(int torindx, Param *paramk, Param *paraml,
tork[2] = delrk[0]*delrj[1] - delrk[1]*delrj[0];
torl[2] = delrj[0]*delrl[1] - delrj[1]*delrl[0];
TT2 = dot3(tork,torl);
rmut = pow((TT2/TT1),2);
rmut = square(TT2/TT1);
if (torindx>=1) {
btt = 1.0 - rmut;
return btt * fc1k * fc1l;
}
else {
btt=paramk->ptork1-TT2/TT1;
btt=paramk->ptork2*pow(btt,2);
btt=paramk->ptork2*square(btt);
return btt * fc1k * fc1l;
}
@ -2930,7 +2932,7 @@ void PairComb3::tor_calc(double r, Param *parami, Param *paramj,
} else {
xtor = -comb_fc(r, parami) * parami->pcross + xcn;
ytor = -comb_fc(r, paramj) * paramj->pcross + ycn;
zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
zcon = 1.0 + square(kconjug) + square(lconjug);
if (xtor < 0.0) xtor = 0.0;
if (ytor < 0.0) ytor = 0.0;
if (zcon < 1.0) zcon = 1.0;
@ -2981,9 +2983,9 @@ void PairComb3::tor_int(int torindx,double xtor, double ytor, double zcon, int l
* pow(ytor,iin3[j][1]) * pow(zcon,iin3[j][2]);
vtor += x;
if (xtor > 1.0e-8 ) dtorx += x*iin3[j][0]/xtor;
if (ytor > 1.0e-8 ) dtory += x*iin3[j][1]/ytor;
if (zcon > 1.0e-8 ) dtorz += x*iin3[j][2]/zcon;
if (xtor > 1.0e-8) dtorx += x*iin3[j][0]/xtor;
if (ytor > 1.0e-8) dtory += x*iin3[j][1]/ytor;
if (zcon > 1.0e-8) dtorz += x*iin3[j][2]/zcon;
}
}
@ -3022,10 +3024,10 @@ void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
fcp1k = comb_fc_d(rik,paramk);
fc1l = comb_fc(rjl,paraml);
fcp1l = comb_fc_d(rjl,paraml);
srmul2 = pow(srmul,2);
srmul2 = square(srmul);
TT1 = rik*rjl*rij*rij*srmu*srmul;
dt1dik = -rmu/pow(srmu,2);
dt1dik = -rmu/square(srmu);
dt1djl = -rmul/srmul2;
tork[0] = delrk[1]*delrj[2] - delrk[2]*delrj[1];
torl[0] = delrj[1]*delrl[2] - delrj[2]*delrl[1];
@ -3051,12 +3053,12 @@ void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
rmut = TT2/TT1;
if (torindx>=1) {
btt = 1.0 - pow(rmut,2);
btt = 1.0 - square(rmut);
AA = -2.0 * ptorr * rmut * fc1k * fc1l / TT1;
}
else {
btt=paramk->ptork1-rmut;
btt=paramk->ptork2*pow(btt,2);
btt=paramk->ptork2*square(btt);
AA = -2.0 * ptorr * paramk->ptork2 *
(paramk->ptork1-rmut) * fc1k * fc1l /TT1;
}
@ -3300,8 +3302,8 @@ double PairComb3::qfo_self(Param *param, double qi)
cmin = cmax = 100.0;
self_d = s1+qi*(2.0*s2+qi*(3.0*s3+qi*4.0*s4));
if (qi < qmin) self_d += 4.0 * cmin * pow((qi-qmin),3);
if (qi > qmax) self_d += 4.0 * cmax * pow((qi-qmax),3);
if (qi < qmin) self_d += 4.0 * cmin * cube((qi-qmin));
if (qi > qmax) self_d += 4.0 * cmax * cube((qi-qmax));
return self_d;
}
@ -3392,10 +3394,10 @@ void PairComb3::qfo_field(Param *parami, Param *paramj, double rsq,
pcmi1 = parami->pcmn1;
pcmi2 = parami->pcmn2;
rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));
rf3i = r3/(square(r3)+cube(pcmi1));
rcf3i = rc3/(square(rc3)+cube(pcmi1));
rf5i = r5/(square(r5)+powint(pcmi2,5));
rcf5i = rc5/(square(rc5)+powint(pcmi2,5));
drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
drcf5i = 5/rc*rcf5i-10*rc4*rcf5i*rcf5i;
@ -3477,7 +3479,7 @@ void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
Di = parami->DU + pow(QUchi,parami->nD); // YYDin
dDi = -parami->nD * parami->bD * pow(QUchi,(parami->nD-1.0)); // YYDiqp
Bsi = parami->aB - pow(QOchi,10); // YYBsin
dBsi = -parami->bB * 10.0 * pow(QOchi,9.0); // YYBsiqp
dBsi = -parami->bB * 10.0 * powint(QOchi,9); // YYBsiqp
}
if (jq < paramj->QL-0.2) {
@ -3492,7 +3494,7 @@ void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
Dj = paramj->DU + pow(QUchj,paramj->nD); // YYDij
dDj = -paramj->nD * paramj->bD * pow(QUchj,(paramj->nD-1.0)); // YYDiqpj
Bsj = paramj->aB - pow(QOchj,10); // YYBsij
dBsj = -paramj->bB * 10.0 * pow(QOchj,9.0); // YYBsiqpj
dBsj = -paramj->bB * 10.0 * powint(QOchj,9); // YYBsiqpj
}
LamDiLamDj = exp(0.5*(parami->lami*Di+paramj->lami*Dj)-rlm1*r);