Compare commits
21 Commits
patch_2Aug
...
coulomb-re
| Author | SHA1 | Date | |
|---|---|---|---|
| a31134a16f | |||
| a89246fbfc | |||
| f739a94eeb | |||
| 9d8d88d020 | |||
| 6410797697 | |||
| 0898ffe336 | |||
| d84ffc5309 | |||
| 124bd99751 | |||
| 85a24bda21 | |||
| d759949d03 | |||
| 25904e80f0 | |||
| cfb2353313 | |||
| bd49f4d665 | |||
| b2ef952187 | |||
| 381d87b79b | |||
| 8778452eaf | |||
| 6ae15d0801 | |||
| 2826449b52 | |||
| c7c6c44add | |||
| 8e145e2b9f | |||
| 0f76bbc718 |
@ -22,21 +22,15 @@
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJClass2CoulLong::PairLJClass2CoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -76,7 +70,7 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double factor_coul,factor_lj;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
@ -130,18 +124,17 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -484,7 +477,7 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r,rinv,r3inv,r6inv,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,forcelj,phicoul,philj;
|
||||
int itable;
|
||||
|
||||
@ -493,18 +486,17 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -145,7 +145,7 @@ void PairBornCoulLongCS::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -145,7 +145,7 @@ void PairBuckCoulLongCS::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -143,7 +143,7 @@ void PairCoulLongCS::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -151,7 +151,7 @@ void PairLJCutCoulLongCS::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -481,7 +481,7 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -22,6 +22,7 @@
|
||||
#include "neigh_list.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
@ -30,16 +31,9 @@
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutDipoleLong::PairLJCutDipoleLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -82,7 +76,7 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag)
|
||||
double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul;
|
||||
double fx,fy,fz,fdx,fdy,fdz,fax,fay,faz;
|
||||
double pdotp,pidotr,pjdotr,pre1,pre2,pre3;
|
||||
double grij,expm2,t,erfc;
|
||||
double grij,expm2,erfc;
|
||||
double g0,g1,g2,b0,b1,b2,b3,d0,d1,d2,d3;
|
||||
double zdix,zdiy,zdiz,zdjx,zdjy,zdjz,zaix,zaiy,zaiz,zajx,zajy,zajz;
|
||||
double g0b1_g1b2_g2b3,g0d1_g1d2_g2d3;
|
||||
@ -112,14 +106,14 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag)
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
pre1 = 2.0 * g_ewald / MY_PIS;
|
||||
pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
|
||||
pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
|
||||
pre2 = 4.0 * g_ewald*g_ewald*g_ewald / MY_PIS;
|
||||
pre3 = 8.0 * g_ewald*g_ewald*g_ewald*g_ewald*g_ewald / MY_PIS;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = atom->q[i];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
@ -140,174 +134,173 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag)
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
rinv = sqrt(r2inv);
|
||||
r2inv = 1.0/rsq;
|
||||
rinv = sqrt(r2inv);
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
if (rsq < cut_coulsq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
|
||||
pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
|
||||
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
|
||||
pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
|
||||
pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
|
||||
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
|
||||
pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
|
||||
|
||||
g0 = qtmp*q[j];
|
||||
g1 = qtmp*pjdotr - q[j]*pidotr + pdotp;
|
||||
g2 = -pidotr*pjdotr;
|
||||
g0 = qtmp*q[j];
|
||||
g1 = qtmp*pjdotr - q[j]*pidotr + pdotp;
|
||||
g2 = -pidotr*pjdotr;
|
||||
|
||||
if (factor_coul > 0.0) {
|
||||
b0 = erfc * rinv;
|
||||
b1 = (b0 + pre1*expm2) * r2inv;
|
||||
b2 = (3.0*b1 + pre2*expm2) * r2inv;
|
||||
b3 = (5.0*b2 + pre3*expm2) * r2inv;
|
||||
if (factor_coul > 0.0) {
|
||||
b0 = erfc * rinv;
|
||||
b1 = (b0 + pre1*expm2) * r2inv;
|
||||
b2 = (3.0*b1 + pre2*expm2) * r2inv;
|
||||
b3 = (5.0*b2 + pre3*expm2) * r2inv;
|
||||
|
||||
g0b1_g1b2_g2b3 = g0*b1 + g1*b2 + g2*b3;
|
||||
fdx = delx * g0b1_g1b2_g2b3 -
|
||||
b1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
|
||||
b2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
|
||||
fdy = dely * g0b1_g1b2_g2b3 -
|
||||
b1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
|
||||
b2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
|
||||
fdz = delz * g0b1_g1b2_g2b3 -
|
||||
b1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
|
||||
b2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
|
||||
g0b1_g1b2_g2b3 = g0*b1 + g1*b2 + g2*b3;
|
||||
fdx = delx * g0b1_g1b2_g2b3 -
|
||||
b1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
|
||||
b2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
|
||||
fdy = dely * g0b1_g1b2_g2b3 -
|
||||
b1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
|
||||
b2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
|
||||
fdz = delz * g0b1_g1b2_g2b3 -
|
||||
b1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
|
||||
b2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
|
||||
|
||||
zdix = delx * (q[j]*b1 + b2*pjdotr) - b1*mu[j][0];
|
||||
zdiy = dely * (q[j]*b1 + b2*pjdotr) - b1*mu[j][1];
|
||||
zdiz = delz * (q[j]*b1 + b2*pjdotr) - b1*mu[j][2];
|
||||
zdjx = delx * (-qtmp*b1 + b2*pidotr) - b1*mu[i][0];
|
||||
zdjy = dely * (-qtmp*b1 + b2*pidotr) - b1*mu[i][1];
|
||||
zdjz = delz * (-qtmp*b1 + b2*pidotr) - b1*mu[i][2];
|
||||
zdix = delx * (q[j]*b1 + b2*pjdotr) - b1*mu[j][0];
|
||||
zdiy = dely * (q[j]*b1 + b2*pjdotr) - b1*mu[j][1];
|
||||
zdiz = delz * (q[j]*b1 + b2*pjdotr) - b1*mu[j][2];
|
||||
zdjx = delx * (-qtmp*b1 + b2*pidotr) - b1*mu[i][0];
|
||||
zdjy = dely * (-qtmp*b1 + b2*pidotr) - b1*mu[i][1];
|
||||
zdjz = delz * (-qtmp*b1 + b2*pidotr) - b1*mu[i][2];
|
||||
|
||||
if (factor_coul < 1.0) {
|
||||
fdx *= factor_coul;
|
||||
fdy *= factor_coul;
|
||||
fdz *= factor_coul;
|
||||
zdix *= factor_coul;
|
||||
zdiy *= factor_coul;
|
||||
zdiz *= factor_coul;
|
||||
zdjx *= factor_coul;
|
||||
zdjy *= factor_coul;
|
||||
zdjz *= factor_coul;
|
||||
}
|
||||
} else {
|
||||
fdx = fdy = fdz = 0.0;
|
||||
zdix = zdiy = zdiz = 0.0;
|
||||
zdjx = zdjy = zdjz = 0.0;
|
||||
}
|
||||
if (factor_coul < 1.0) {
|
||||
fdx *= factor_coul;
|
||||
fdy *= factor_coul;
|
||||
fdz *= factor_coul;
|
||||
zdix *= factor_coul;
|
||||
zdiy *= factor_coul;
|
||||
zdiz *= factor_coul;
|
||||
zdjx *= factor_coul;
|
||||
zdjy *= factor_coul;
|
||||
zdjz *= factor_coul;
|
||||
}
|
||||
} else {
|
||||
fdx = fdy = fdz = 0.0;
|
||||
zdix = zdiy = zdiz = 0.0;
|
||||
zdjx = zdjy = zdjz = 0.0;
|
||||
}
|
||||
|
||||
if (factor_coul < 1.0) {
|
||||
d0 = (erfc - 1.0) * rinv;
|
||||
d1 = (d0 + pre1*expm2) * r2inv;
|
||||
d2 = (3.0*d1 + pre2*expm2) * r2inv;
|
||||
d3 = (5.0*d2 + pre3*expm2) * r2inv;
|
||||
if (factor_coul < 1.0) {
|
||||
d0 = (erfc - 1.0) * rinv;
|
||||
d1 = (d0 + pre1*expm2) * r2inv;
|
||||
d2 = (3.0*d1 + pre2*expm2) * r2inv;
|
||||
d3 = (5.0*d2 + pre3*expm2) * r2inv;
|
||||
|
||||
g0d1_g1d2_g2d3 = g0*d1 + g1*d2 + g2*d3;
|
||||
fax = delx * g0d1_g1d2_g2d3 -
|
||||
d1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
|
||||
d2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
|
||||
fay = dely * g0d1_g1d2_g2d3 -
|
||||
d1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
|
||||
d2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
|
||||
faz = delz * g0d1_g1d2_g2d3 -
|
||||
d1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
|
||||
d2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
|
||||
g0d1_g1d2_g2d3 = g0*d1 + g1*d2 + g2*d3;
|
||||
fax = delx * g0d1_g1d2_g2d3 -
|
||||
d1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
|
||||
d2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
|
||||
fay = dely * g0d1_g1d2_g2d3 -
|
||||
d1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
|
||||
d2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
|
||||
faz = delz * g0d1_g1d2_g2d3 -
|
||||
d1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
|
||||
d2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
|
||||
|
||||
zaix = delx * (q[j]*d1 + d2*pjdotr) - d1*mu[j][0];
|
||||
zaiy = dely * (q[j]*d1 + d2*pjdotr) - d1*mu[j][1];
|
||||
zaiz = delz * (q[j]*d1 + d2*pjdotr) - d1*mu[j][2];
|
||||
zajx = delx * (-qtmp*d1 + d2*pidotr) - d1*mu[i][0];
|
||||
zajy = dely * (-qtmp*d1 + d2*pidotr) - d1*mu[i][1];
|
||||
zajz = delz * (-qtmp*d1 + d2*pidotr) - d1*mu[i][2];
|
||||
zaix = delx * (q[j]*d1 + d2*pjdotr) - d1*mu[j][0];
|
||||
zaiy = dely * (q[j]*d1 + d2*pjdotr) - d1*mu[j][1];
|
||||
zaiz = delz * (q[j]*d1 + d2*pjdotr) - d1*mu[j][2];
|
||||
zajx = delx * (-qtmp*d1 + d2*pidotr) - d1*mu[i][0];
|
||||
zajy = dely * (-qtmp*d1 + d2*pidotr) - d1*mu[i][1];
|
||||
zajz = delz * (-qtmp*d1 + d2*pidotr) - d1*mu[i][2];
|
||||
|
||||
if (factor_coul > 0.0) {
|
||||
facm1 = 1.0 - factor_coul;
|
||||
fax *= facm1;
|
||||
fay *= facm1;
|
||||
faz *= facm1;
|
||||
zaix *= facm1;
|
||||
zaiy *= facm1;
|
||||
zaiz *= facm1;
|
||||
zajx *= facm1;
|
||||
zajy *= facm1;
|
||||
zajz *= facm1;
|
||||
}
|
||||
} else {
|
||||
fax = fay = faz = 0.0;
|
||||
zaix = zaiy = zaiz = 0.0;
|
||||
zajx = zajy = zajz = 0.0;
|
||||
}
|
||||
if (factor_coul > 0.0) {
|
||||
facm1 = 1.0 - factor_coul;
|
||||
fax *= facm1;
|
||||
fay *= facm1;
|
||||
faz *= facm1;
|
||||
zaix *= facm1;
|
||||
zaiy *= facm1;
|
||||
zaiz *= facm1;
|
||||
zajx *= facm1;
|
||||
zajy *= facm1;
|
||||
zajz *= facm1;
|
||||
}
|
||||
} else {
|
||||
fax = fay = faz = 0.0;
|
||||
zaix = zaiy = zaiz = 0.0;
|
||||
zajx = zajy = zajz = 0.0;
|
||||
}
|
||||
|
||||
forcecoulx = fdx + fax;
|
||||
forcecouly = fdy + fay;
|
||||
forcecoulz = fdz + faz;
|
||||
forcecoulx = fdx + fax;
|
||||
forcecouly = fdy + fay;
|
||||
forcecoulz = fdz + faz;
|
||||
|
||||
tixcoul = mu[i][1]*(zdiz + zaiz) - mu[i][2]*(zdiy + zaiy);
|
||||
tiycoul = mu[i][2]*(zdix + zaix) - mu[i][0]*(zdiz + zaiz);
|
||||
tizcoul = mu[i][0]*(zdiy + zaiy) - mu[i][1]*(zdix + zaix);
|
||||
tjxcoul = mu[j][1]*(zdjz + zajz) - mu[j][2]*(zdjy + zajy);
|
||||
tjycoul = mu[j][2]*(zdjx + zajx) - mu[j][0]*(zdjz + zajz);
|
||||
tjzcoul = mu[j][0]*(zdjy + zajy) - mu[j][1]*(zdjx + zajx);
|
||||
tixcoul = mu[i][1]*(zdiz + zaiz) - mu[i][2]*(zdiy + zaiy);
|
||||
tiycoul = mu[i][2]*(zdix + zaix) - mu[i][0]*(zdiz + zaiz);
|
||||
tizcoul = mu[i][0]*(zdiy + zaiy) - mu[i][1]*(zdix + zaix);
|
||||
tjxcoul = mu[j][1]*(zdjz + zajz) - mu[j][2]*(zdjy + zajy);
|
||||
tjycoul = mu[j][2]*(zdjx + zajx) - mu[j][0]*(zdjz + zajz);
|
||||
tjzcoul = mu[j][0]*(zdjy + zajy) - mu[j][1]*(zdjx + zajx);
|
||||
|
||||
} else {
|
||||
forcecoulx = forcecouly = forcecoulz = 0.0;
|
||||
tixcoul = tiycoul = tizcoul = 0.0;
|
||||
tjxcoul = tjycoul = tjzcoul = 0.0;
|
||||
}
|
||||
|
||||
// LJ interaction
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
fforce = factor_lj * forcelj*r2inv;
|
||||
} else fforce = 0.0;
|
||||
|
||||
// total force
|
||||
|
||||
fx = qqrd2e*forcecoulx + delx*fforce;
|
||||
fy = qqrd2e*forcecouly + dely*fforce;
|
||||
fz = qqrd2e*forcecoulz + delz*fforce;
|
||||
|
||||
// force & torque accumulation
|
||||
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
torque[i][0] += qqrd2e*tixcoul;
|
||||
torque[i][1] += qqrd2e*tiycoul;
|
||||
torque[i][2] += qqrd2e*tizcoul;
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= fx;
|
||||
f[j][1] -= fy;
|
||||
f[j][2] -= fz;
|
||||
torque[j][0] += qqrd2e*tjxcoul;
|
||||
torque[j][1] += qqrd2e*tjycoul;
|
||||
torque[j][2] += qqrd2e*tjzcoul;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq && factor_coul > 0.0) {
|
||||
ecoul = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
|
||||
if (factor_coul < 1.0) {
|
||||
ecoul *= factor_coul;
|
||||
ecoul += (1-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
|
||||
} else {
|
||||
forcecoulx = forcecouly = forcecoulz = 0.0;
|
||||
tixcoul = tiycoul = tizcoul = 0.0;
|
||||
tjxcoul = tjycoul = tjzcoul = 0.0;
|
||||
}
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
// LJ interaction
|
||||
|
||||
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
|
||||
evdwl,ecoul,fx,fy,fz,delx,dely,delz);
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
fforce = factor_lj * forcelj*r2inv;
|
||||
} else fforce = 0.0;
|
||||
|
||||
// total force
|
||||
|
||||
fx = qqrd2e*forcecoulx + delx*fforce;
|
||||
fy = qqrd2e*forcecouly + dely*fforce;
|
||||
fz = qqrd2e*forcecoulz + delz*fforce;
|
||||
|
||||
// force & torque accumulation
|
||||
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
torque[i][0] += qqrd2e*tixcoul;
|
||||
torque[i][1] += qqrd2e*tiycoul;
|
||||
torque[i][2] += qqrd2e*tizcoul;
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= fx;
|
||||
f[j][1] -= fy;
|
||||
f[j][2] -= fz;
|
||||
torque[j][0] += qqrd2e*tjxcoul;
|
||||
torque[j][1] += qqrd2e*tjycoul;
|
||||
torque[j][2] += qqrd2e*tjzcoul;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq && factor_coul > 0.0) {
|
||||
ecoul = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
|
||||
if (factor_coul < 1.0) {
|
||||
ecoul *= factor_coul;
|
||||
ecoul += (1-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
|
||||
}
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
|
||||
evdwl,ecoul,fx,fy,fz,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -260,7 +260,7 @@ void PairCoulLongGPU::cpu_compute(int start, int inum, int eflag,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -282,7 +282,7 @@ void PairLJCharmmCoulLongGPU::cpu_compute(int start, int inum, int eflag,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -278,7 +278,7 @@ void PairLJCutCoulLongGPU::cpu_compute(int start, int inum, int eflag,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -250,7 +250,7 @@ void PairLJCutCoulMSMGPU::cpu_compute(int start, int inum, int eflag, int vflag,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -279,7 +279,7 @@ void PairLJSDKCoulLongGPU::cpu_compute(int start, int inum, int *ilist,
|
||||
rsq_lookup.f = rsq;
|
||||
int itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
const double fraction = (rsq_lookup.f - rtable[itable]) *
|
||||
const double fraction = (rsq - rtable[itable]) *
|
||||
drtable[itable];
|
||||
const double table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
|
||||
@ -315,7 +315,7 @@ void PairTableGPU::cpu_compute(int start, int inum, int eflag, int vflag,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & tb->nmask;
|
||||
itable >>= tb->nshiftbits;
|
||||
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
|
||||
fraction = (rsq - tb->rsq[itable]) * tb->drsq[itable];
|
||||
value = tb->f[itable] + fraction*tb->df[itable];
|
||||
fpair = factor_lj * value;
|
||||
}
|
||||
|
||||
@ -30,26 +30,20 @@
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "atom_masks.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
@ -236,7 +230,7 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
|
||||
F_FLOAT forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -248,12 +242,11 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT rinv = 1.0/r;
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
|
||||
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
F_FLOAT forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
return forcecoul*rinv*rinv;
|
||||
@ -274,7 +267,7 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
|
||||
F_FLOAT ecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -286,9 +279,8 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
F_FLOAT ecoul = prefactor * erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -30,26 +30,20 @@
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "atom_masks.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
@ -200,7 +194,7 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
|
||||
F_FLOAT forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -212,12 +206,11 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT rinv = 1.0/r;
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
|
||||
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
F_FLOAT forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
return forcecoul*rinv*rinv;
|
||||
@ -238,7 +231,7 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
|
||||
F_FLOAT ecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -250,9 +243,8 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
F_FLOAT ecoul = prefactor * erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -41,15 +41,6 @@ using namespace MathConst;
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
|
||||
@ -41,15 +41,6 @@ using namespace MathConst;
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
|
||||
@ -30,26 +30,20 @@
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "atom_masks.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
@ -265,7 +259,7 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
|
||||
F_FLOAT forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -277,12 +271,11 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT rinv = 1.0/r;
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
|
||||
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
F_FLOAT forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
return forcecoul*rinv*rinv;
|
||||
@ -302,7 +295,7 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
|
||||
F_FLOAT ecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -314,9 +307,8 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
F_FLOAT ecoul = prefactor * erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -26,26 +26,20 @@
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "atom_masks.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
@ -216,7 +210,7 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
|
||||
F_FLOAT forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -228,12 +222,11 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT rinv = 1.0/r;
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
|
||||
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
F_FLOAT forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
return forcecoul*rinv*rinv;
|
||||
@ -274,7 +267,7 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
|
||||
F_FLOAT ecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -286,9 +279,8 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
F_FLOAT ecoul = prefactor * erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -26,26 +26,20 @@
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "atom_masks.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define KOKKOS_CUDA_MAX_THREADS 256
|
||||
#define KOKKOS_CUDA_MIN_BLOCKS 8
|
||||
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType>
|
||||
@ -215,7 +209,7 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
|
||||
F_FLOAT forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -227,12 +221,11 @@ compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT rinv = 1.0/r;
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
|
||||
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
F_FLOAT forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
return forcecoul*rinv*rinv;
|
||||
@ -271,7 +264,7 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT fraction = (rsq - d_rtable[itable]) * d_drtable[itable];
|
||||
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
|
||||
F_FLOAT ecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -283,9 +276,8 @@ compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
|
||||
} else {
|
||||
const F_FLOAT r = sqrt(rsq);
|
||||
const F_FLOAT grij = g_ewald * r;
|
||||
const F_FLOAT expm2 = exp(-grij*grij);
|
||||
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const F_FLOAT expm2 = expmsq(grij);
|
||||
const F_FLOAT erfc = my_erfcx(grij) * expm2;
|
||||
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
F_FLOAT ecoul = prefactor * erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -220,7 +220,7 @@ compute_fpair(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, c
|
||||
rsq_lookup.f = rsq;
|
||||
int itable = rsq_lookup.i & d_table_const.nmask(tidx);
|
||||
itable >>= d_table_const.nshiftbits(tidx);
|
||||
const double fraction = (rsq_lookup.f - d_table_const.rsq(tidx,itable)) * d_table_const.drsq(tidx,itable);
|
||||
const double fraction = (rsq - d_table_const.rsq(tidx,itable)) * d_table_const.drsq(tidx,itable);
|
||||
fpair = d_table_const.f(tidx,itable) + fraction*d_table_const.df(tidx,itable);
|
||||
}
|
||||
return fpair;
|
||||
@ -254,7 +254,7 @@ compute_evdwl(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, c
|
||||
rsq_lookup.f = rsq;
|
||||
int itable = rsq_lookup.i & d_table_const.nmask(tidx);
|
||||
itable >>= d_table_const.nshiftbits(tidx);
|
||||
const double fraction = (rsq_lookup.f - d_table_const.rsq(tidx,itable)) * d_table_const.drsq(tidx,itable);
|
||||
const double fraction = (rsq - d_table_const.rsq(tidx,itable)) * d_table_const.drsq(tidx,itable);
|
||||
evdwl = d_table_const.e(tidx,itable) + fraction*d_table_const.de(tidx,itable);
|
||||
}
|
||||
return evdwl;
|
||||
|
||||
@ -26,21 +26,15 @@
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBornCoulLong::PairBornCoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -82,7 +76,7 @@ void PairBornCoulLong::compute(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double r,rexp;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
@ -136,18 +130,17 @@ void PairBornCoulLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -514,7 +507,7 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r6inv,r,rexp,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,forceborn,phicoul,phiborn;
|
||||
int itable;
|
||||
|
||||
@ -523,18 +516,17 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -22,21 +22,15 @@
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBuckCoulLong::PairBuckCoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -77,7 +71,7 @@ void PairBuckCoulLong::compute(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double r,rexp;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
@ -126,22 +120,22 @@ void PairBuckCoulLong::compute(int eflag, int vflag)
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -486,7 +480,7 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r6inv,r,rexp,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,forcebuck,phicoul,phibuck;
|
||||
int itable;
|
||||
|
||||
@ -495,18 +489,17 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -24,23 +24,19 @@
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -74,7 +70,7 @@ void PairCoulLong::compute(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,forcecoul,factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
|
||||
@ -124,18 +120,17 @@ void PairCoulLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -343,7 +338,7 @@ double PairCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,phicoul;
|
||||
int itable;
|
||||
|
||||
@ -351,18 +346,17 @@ double PairCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -121,7 +121,7 @@ void PairCoulMSM::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -192,7 +192,7 @@ double PairCoulMSM::single(int i, int j, int itype, int jtype,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -30,18 +30,14 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -90,7 +86,7 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double philj,switch1,switch2;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
@ -144,18 +140,17 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -397,7 +392,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double philj,switch1,switch2;
|
||||
double rsw;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
@ -453,11 +448,10 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2 - 1.0);
|
||||
if (rsq > cut_in_off_sq) {
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (r - cut_in_off)/cut_in_diff;
|
||||
@ -476,7 +470,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -546,7 +540,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
|
||||
if (vflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
table = vtable[itable] + fraction*dvtable[itable];
|
||||
@ -940,7 +934,7 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r6inv,r,grij,expm2,erfc,prefactor;
|
||||
double switch1,switch2,fraction,table,forcecoul,forcelj,phicoul,philj;
|
||||
int itable;
|
||||
|
||||
@ -949,18 +943,17 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -144,7 +144,7 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -342,7 +342,7 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -484,7 +484,7 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -30,21 +30,15 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulLong::PairLJCutCoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -85,7 +79,7 @@ void PairLJCutCoulLong::compute(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
|
||||
@ -139,18 +133,17 @@ void PairLJCutCoulLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -391,7 +384,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double rsw;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
@ -453,11 +446,10 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2 - 1.0);
|
||||
if (rsq > cut_in_off_sq) {
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (r - cut_in_off)/cut_in_diff;
|
||||
@ -476,7 +468,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -534,7 +526,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
|
||||
if (vflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
table = vtable[itable] + fraction*dvtable[itable];
|
||||
@ -906,7 +898,7 @@ double PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r6inv,r,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,forcelj,phicoul,philj;
|
||||
int itable;
|
||||
|
||||
@ -915,18 +907,17 @@ double PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup_single;
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -145,7 +145,7 @@ void PairLJCutCoulMSM::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -329,7 +329,7 @@ void PairLJCutCoulMSM::compute_outer(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -442,7 +442,7 @@ double PairLJCutCoulMSM::single(int i, int j, int itype, int jtype,
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -33,18 +33,14 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -87,7 +83,7 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,cforce;
|
||||
double factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double fO[3],fH[3],fd[3],v[6];
|
||||
double *x1,*x2,*xH1,*xH2;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
@ -260,11 +256,10 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
@ -273,7 +268,7 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -309,7 +309,7 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -33,18 +33,14 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -83,7 +79,7 @@ void PairTIP4PLong::compute(int eflag, int vflag)
|
||||
double fraction,table;
|
||||
double r,r2inv,forcecoul,cforce;
|
||||
double factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double fO[3],fH[3],fd[3],v[6];
|
||||
double *x1,*x2,*xH1,*xH2;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
@ -228,11 +224,10 @@ void PairTIP4PLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
@ -241,7 +236,7 @@ void PairTIP4PLong::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -1209,23 +1209,24 @@ double PPPM::compute_qopt()
|
||||
sum2 = 0.0;
|
||||
sum3 = 0.0;
|
||||
sum4 = 0.0;
|
||||
const double inv_gew = 1.0/g_ewald;
|
||||
for (nx = -2; nx <= 2; nx++) {
|
||||
qx = unitkx*(kper+nx_pppm*nx);
|
||||
sx = exp(-0.25*square(qx/g_ewald));
|
||||
sx = expmsq(0.5*qx*inv_gew);
|
||||
argx = 0.5*qx*xprd/nx_pppm;
|
||||
wx = powsinxx(argx,twoorder);
|
||||
qx *= qx;
|
||||
|
||||
for (ny = -2; ny <= 2; ny++) {
|
||||
qy = unitky*(lper+ny_pppm*ny);
|
||||
sy = exp(-0.25*square(qy/g_ewald));
|
||||
sy = expmsq(0.5*qy*inv_gew);
|
||||
argy = 0.5*qy*yprd/ny_pppm;
|
||||
wy = powsinxx(argy,twoorder);
|
||||
qy *= qy;
|
||||
|
||||
for (nz = -2; nz <= 2; nz++) {
|
||||
qz = unitkz*(mper+nz_pppm*nz);
|
||||
sz = exp(-0.25*square(qz/g_ewald));
|
||||
sz = expmsq(0.5*qz*inv_gew);
|
||||
argz = 0.5*qz*zprd_slab/nz_pppm;
|
||||
wz = powsinxx(argz,twoorder);
|
||||
qz *= qz;
|
||||
@ -1297,7 +1298,7 @@ double PPPM::newton_raphson_f()
|
||||
double zprd = domain->zprd;
|
||||
bigint natoms = atom->natoms;
|
||||
|
||||
double df_rspace = 2.0*q2*exp(-g_ewald*g_ewald*cutoff*cutoff) /
|
||||
double df_rspace = 2.0*q2*expmsq(g_ewald*cutoff) /
|
||||
sqrt(natoms*cutoff*xprd*yprd*zprd);
|
||||
|
||||
double df_kspace = compute_df_kspace();
|
||||
@ -1338,7 +1339,7 @@ double PPPM::final_accuracy()
|
||||
|
||||
double df_kspace = compute_df_kspace();
|
||||
double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd);
|
||||
double df_rspace = 2.0 * q2_over_sqrt * exp(-g_ewald*g_ewald*cutoff*cutoff);
|
||||
double df_rspace = 2.0 * q2_over_sqrt * expmsq(g_ewald*cutoff);
|
||||
double df_table = estimate_table_accuracy(q2_over_sqrt,df_rspace);
|
||||
double estimated_accuracy = sqrt(df_kspace*df_kspace + df_rspace*df_rspace +
|
||||
df_table*df_table);
|
||||
@ -1579,22 +1580,23 @@ void PPPM::compute_gf_ik()
|
||||
numerator = 12.5663706/sqk;
|
||||
denominator = gf_denom(snx,sny,snz);
|
||||
sum1 = 0.0;
|
||||
const double inv_gew = 1.0/g_ewald;
|
||||
|
||||
for (nx = -nbx; nx <= nbx; nx++) {
|
||||
qx = unitkx*(kper+nx_pppm*nx);
|
||||
sx = exp(-0.25*square(qx/g_ewald));
|
||||
sx = expmsq(0.5*qx*inv_gew);
|
||||
argx = 0.5*qx*xprd/nx_pppm;
|
||||
wx = powsinxx(argx,twoorder);
|
||||
|
||||
for (ny = -nby; ny <= nby; ny++) {
|
||||
qy = unitky*(lper+ny_pppm*ny);
|
||||
sy = exp(-0.25*square(qy/g_ewald));
|
||||
sy = expmsq(0.5*qy*inv_gew);
|
||||
argy = 0.5*qy*yprd/ny_pppm;
|
||||
wy = powsinxx(argy,twoorder);
|
||||
|
||||
for (nz = -nbz; nz <= nbz; nz++) {
|
||||
qz = unitkz*(mper+nz_pppm*nz);
|
||||
sz = exp(-0.25*square(qz/g_ewald));
|
||||
sz = expmsq(0.5*qz*inv_gew);
|
||||
argz = 0.5*qz*zprd_slab/nz_pppm;
|
||||
wz = powsinxx(argz,twoorder);
|
||||
|
||||
@ -1662,6 +1664,7 @@ void PPPM::compute_gf_ik_triclinic()
|
||||
numerator = 12.5663706/sqk;
|
||||
denominator = gf_denom(snx,sny,snz);
|
||||
sum1 = 0.0;
|
||||
const double inv_gew = 1.0/g_ewald;
|
||||
|
||||
for (nx = -nbx; nx <= nbx; nx++) {
|
||||
argx = MY_PI*kper/nx_pppm + MY_PI*nx;
|
||||
@ -1682,13 +1685,13 @@ void PPPM::compute_gf_ik_triclinic()
|
||||
x2lamdaT(&b[0],&b[0]);
|
||||
|
||||
qx = unitk_lamda[0]+b[0];
|
||||
sx = exp(-0.25*square(qx/g_ewald));
|
||||
sx = expmsq(0.5*qx*inv_gew);
|
||||
|
||||
qy = unitk_lamda[1]+b[1];
|
||||
sy = exp(-0.25*square(qy/g_ewald));
|
||||
sy = expmsq(0.5*qy*inv_gew);
|
||||
|
||||
qz = unitk_lamda[2]+b[2];
|
||||
sz = exp(-0.25*square(qz/g_ewald));
|
||||
sz = expmsq(0.5*qz*inv_gew);
|
||||
|
||||
dot1 = unitk_lamda[0]*qx + unitk_lamda[1]*qy + unitk_lamda[2]*qz;
|
||||
dot2 = qx*qx+qy*qy+qz*qz;
|
||||
@ -1725,6 +1728,7 @@ void PPPM::compute_gf_ad()
|
||||
int k,l,m,n,kper,lper,mper;
|
||||
|
||||
const int twoorder = 2*order;
|
||||
const double inv_gew = 1.0/g_ewald;
|
||||
|
||||
for (int i = 0; i < 6; i++) sf_coeff[i] = 0.0;
|
||||
|
||||
@ -1733,7 +1737,7 @@ void PPPM::compute_gf_ad()
|
||||
mper = m - nz_pppm*(2*m/nz_pppm);
|
||||
qz = unitkz*mper;
|
||||
snz = square(sin(0.5*qz*zprd_slab/nz_pppm));
|
||||
sz = exp(-0.25*square(qz/g_ewald));
|
||||
sz = expmsq(0.5*qz*inv_gew);
|
||||
argz = 0.5*qz*zprd_slab/nz_pppm;
|
||||
wz = powsinxx(argz,twoorder);
|
||||
|
||||
@ -1741,7 +1745,7 @@ void PPPM::compute_gf_ad()
|
||||
lper = l - ny_pppm*(2*l/ny_pppm);
|
||||
qy = unitky*lper;
|
||||
sny = square(sin(0.5*qy*yprd/ny_pppm));
|
||||
sy = exp(-0.25*square(qy/g_ewald));
|
||||
sy = expmsq(0.5*qy*inv_gew);
|
||||
argy = 0.5*qy*yprd/ny_pppm;
|
||||
wy = powsinxx(argy,twoorder);
|
||||
|
||||
@ -1749,7 +1753,7 @@ void PPPM::compute_gf_ad()
|
||||
kper = k - nx_pppm*(2*k/nx_pppm);
|
||||
qx = unitkx*kper;
|
||||
snx = square(sin(0.5*qx*xprd/nx_pppm));
|
||||
sx = exp(-0.25*square(qx/g_ewald));
|
||||
sx = expmsq(0.5*qx*inv_gew);
|
||||
argx = 0.5*qx*xprd/nx_pppm;
|
||||
wx = powsinxx(argx,twoorder);
|
||||
|
||||
|
||||
@ -29,21 +29,15 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairNMCutCoulLong::PairNMCutCoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -51,6 +45,7 @@ PairNMCutCoulLong::PairNMCutCoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
ewaldflag = pppmflag = 1;
|
||||
ftable = NULL;
|
||||
qdist = 0.0;
|
||||
writedata = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -85,7 +80,7 @@ void PairNMCutCoulLong::compute(int eflag, int vflag)
|
||||
double fraction,table;
|
||||
double r,r2inv,factor_coul,factor_lj;
|
||||
double forcecoul,forcenm,rminv,rninv;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
|
||||
@ -139,18 +134,17 @@ void PairNMCutCoulLong::compute(int eflag, int vflag)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -512,7 +506,7 @@ double PairNMCutCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,forcenm,phicoul,phinm;
|
||||
int itable;
|
||||
|
||||
@ -521,18 +515,17 @@ double PairNMCutCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup_single;
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -23,17 +23,13 @@
|
||||
#include "pair_lj_charmm_coul_long_opt.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define EWALD_A1 0.254829592
|
||||
#define EWALD_A2 -0.284496736
|
||||
#define EWALD_A3 1.421413741
|
||||
#define EWALD_A4 -1.453152027
|
||||
#define EWALD_A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -76,7 +72,7 @@ void PairLJCharmmCoulLongOpt::eval()
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype,itable,sbindex;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double philj,switch1,switch2;
|
||||
|
||||
double rsq;
|
||||
@ -156,19 +152,16 @@ void PairLJCharmmCoulLongOpt::eval()
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t *
|
||||
(EWALD_A1+t*(EWALD_A2+t*(EWALD_A3+t*(EWALD_A4+t*EWALD_A5)))) *
|
||||
expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * tmp_coef3/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = tmp_coef3 * table;
|
||||
}
|
||||
@ -245,22 +238,17 @@ void PairLJCharmmCoulLongOpt::eval()
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t *
|
||||
(EWALD_A1+t*(EWALD_A2+t*(EWALD_A3+t*(EWALD_A4+t*EWALD_A5)))) *
|
||||
expm2;
|
||||
prefactor = qqrd2e * tmp_coef3/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = tmp_coef3 * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -15,17 +15,13 @@
|
||||
#include "pair_lj_cut_coul_long_opt.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -79,7 +75,7 @@ void PairLJCutCoulLongOpt::eval()
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
|
||||
@ -132,18 +128,17 @@ void PairLJCutCoulLongOpt::eval()
|
||||
if (!CTABLE || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -22,18 +22,14 @@
|
||||
#include "force.h"
|
||||
#include "error.h"
|
||||
#include "memory.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -110,7 +106,7 @@ void PairLJCutTIP4PLongOpt::eval()
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
|
||||
double factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double v[6];
|
||||
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
|
||||
const double *x1,*x2,*xH1,*xH2;
|
||||
@ -272,11 +268,10 @@ void PairLJCutTIP4PLongOpt::eval()
|
||||
if (CTABLE || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
@ -285,7 +280,7 @@ void PairLJCutTIP4PLongOpt::eval()
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -28,8 +28,7 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
@ -37,17 +36,10 @@
|
||||
#include "lj_sdk_common.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
using namespace LJSDKParms;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJSDKCoulLong::PairLJSDKCoulLong(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -118,7 +110,7 @@ void PairLJSDKCoulLong::eval()
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = atom->f;
|
||||
@ -172,11 +164,10 @@ void PairLJSDKCoulLong::eval()
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (EFLAG) ecoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
@ -187,7 +178,7 @@ void PairLJSDKCoulLong::eval()
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] *
|
||||
@ -557,7 +548,7 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,grij,expm2,t,erfc,prefactor;
|
||||
double r2inv,r,grij,expm2,erfc,prefactor;
|
||||
double fraction,table,forcecoul,forcelj,phicoul,philj;
|
||||
int itable;
|
||||
|
||||
@ -568,11 +559,10 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
phicoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
@ -583,7 +573,7 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
|
||||
@ -155,7 +155,7 @@ void PairLJSDKCoulMSM::eval_msm()
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] *
|
||||
@ -253,7 +253,7 @@ double PairLJSDKCoulMSM::single(int i, int j, int itype, int jtype,
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
|
||||
@ -190,7 +190,7 @@ void PairTableRX::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & tb->nmask;
|
||||
itable >>= tb->nshiftbits;
|
||||
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
|
||||
fraction = (rsq - tb->rsq[itable]) * tb->drsq[itable];
|
||||
value = tb->f[itable] + fraction*tb->df[itable];
|
||||
fpair = factor_lj * value;
|
||||
}
|
||||
@ -1107,7 +1107,7 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & tb->nmask;
|
||||
itable >>= tb->nshiftbits;
|
||||
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
|
||||
fraction = (rsq - tb->rsq[itable]) * tb->drsq[itable];
|
||||
value = tb->f[itable] + fraction*tb->df[itable];
|
||||
fforce = factor_lj * value;
|
||||
}
|
||||
|
||||
@ -170,7 +170,7 @@ void PairLJCutTholeLong::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qi*qj * table;
|
||||
if (factor_coul < 1.0) {
|
||||
@ -625,7 +625,7 @@ double PairLJCutTholeLong::single(int i, int j, int itype, int jtype,
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -319,7 +319,7 @@ void PairBuckCoulLongIntel::eval(const int offload, const int vflag,
|
||||
float rsq_lookup = rsq;
|
||||
const int itable = (__intel_castf32_u32(rsq_lookup) &
|
||||
ncoulmask) >> ncoulshiftbits;
|
||||
const flt_t fraction = (rsq_lookup - table[itable].r) *
|
||||
const flt_t fraction = (rsq - table[itable].r) *
|
||||
table[itable].dr;
|
||||
|
||||
const flt_t tablet = table[itable].f +
|
||||
|
||||
@ -354,10 +354,9 @@ void PairLJCharmmCoulLongIntel::eval(const int offload, const int vflag,
|
||||
|
||||
#ifdef INTEL_ALLOW_TABLE
|
||||
} else {
|
||||
float rsq_lookup = rsq;
|
||||
const int itable = (__intel_castf32_u32(rsq_lookup) &
|
||||
ncoulmask) >> ncoulshiftbits;
|
||||
const flt_t fraction = (rsq_lookup - table[itable].r) *
|
||||
const flt_t fraction = (rsq - table[itable].r) *
|
||||
table[itable].dr;
|
||||
|
||||
const flt_t tablet = table[itable].f +
|
||||
|
||||
@ -316,7 +316,7 @@ void PairLJCutCoulLongIntel::eval(const int offload, const int vflag,
|
||||
float rsq_lookup = rsq;
|
||||
const int itable = (__intel_castf32_u32(rsq_lookup) &
|
||||
ncoulmask) >> ncoulshiftbits;
|
||||
const flt_t fraction = (rsq_lookup - table[itable].r) *
|
||||
const flt_t fraction = (rsq - table[itable].r) *
|
||||
table[itable].dr;
|
||||
|
||||
const flt_t tablet = table[itable].f +
|
||||
|
||||
@ -17,19 +17,15 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -89,8 +85,8 @@ void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rsq,r2inv,r6inv,r,rexp,forcecoul,forceborn,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double grij,expm2,prefactor,erfc,table,fraction;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh,itable;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
@ -136,19 +132,34 @@ void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r = sqrt(rsq);
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);
|
||||
forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
|
||||
@ -168,7 +179,12 @@ void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
if (EFLAG) {
|
||||
if (rsq < cut_coulsq) {
|
||||
ecoul = prefactor*erfc;
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j] * table;
|
||||
}
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
|
||||
@ -17,19 +17,15 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -90,8 +86,8 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rsq,r2inv,r6inv,r,rexp,forcecoul,forcebuck,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double grij,expm2,prefactor,erfc,fraction,table;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh,itable;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
@ -137,19 +133,34 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r = sqrt(rsq);
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
rexp = exp(-r*rhoinv[itype][jtype]);
|
||||
forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
|
||||
@ -168,7 +179,12 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
if (EFLAG) {
|
||||
if (rsq < cut_coulsq) {
|
||||
ecoul = prefactor*erfc;
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j] * table;
|
||||
}
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
|
||||
@ -17,22 +17,17 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
#include "math_const.h"
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulDSFOMP::PairCoulDSFOMP(LAMMPS *lmp) :
|
||||
@ -91,7 +86,7 @@ void PairCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
int i,j,ii,jj,jnum;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double r,rsq,r2inv,forcecoul,factor_coul;
|
||||
double prefactor,erfcc,erfcd,t;
|
||||
double prefactor,erfcc,erfcd,grij;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
@ -141,9 +136,9 @@ void PairCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e*qtmp*q[j]/r;
|
||||
erfcd = exp(-alpha*alpha*rsq);
|
||||
t = 1.0 / (1.0 + EWALD_P*alpha*r);
|
||||
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
|
||||
grij = alpha*r;
|
||||
erfcd = expmsq(grij);
|
||||
erfcc = my_erfcx(grij) * erfcd;
|
||||
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
|
||||
r*f_shift) * r;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -17,19 +17,15 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -91,7 +87,7 @@ void PairCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,rsq,forcecoul,factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
@ -139,18 +135,17 @@ void PairCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -145,7 +145,7 @@ void PairCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -17,11 +17,15 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -132,21 +136,12 @@ void PairLJCharmmCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
const double A1 = 0.254829592;
|
||||
const double A2 = -0.284496736;
|
||||
const double A3 = 1.421413741;
|
||||
const double A4 = -1.453152027;
|
||||
const double A5 = 1.061405429;
|
||||
const double EWALD_F = 1.12837917;
|
||||
const double INV_EWALD_P = 1.0/0.3275911;
|
||||
|
||||
const double r = sqrt(rsq);
|
||||
const double grij = g_ewald * r;
|
||||
const double expm2 = exp(-grij*grij);
|
||||
const double t = INV_EWALD_P / (INV_EWALD_P + grij);
|
||||
const double erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const double expm2 = expmsq(grij);
|
||||
const double erfc = my_erfcx(grij) * expm2;
|
||||
const double prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (EFLAG) ecoul = prefactor*erfc;
|
||||
if (sbindex) {
|
||||
const double adjust = (1.0-special_coul[sbindex])*prefactor;
|
||||
@ -157,7 +152,7 @@ void PairLJCharmmCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
const double fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
const double table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
|
||||
|
||||
@ -150,7 +150,7 @@ void PairLJCharmmCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
const double fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
const double table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
|
||||
|
||||
@ -17,19 +17,15 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -89,7 +85,7 @@ void PairLJClass2CoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -140,11 +136,10 @@ void PairLJClass2CoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (rsq < cut_coulsq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
|
||||
@ -17,22 +17,17 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
#include "math_const.h"
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulDSFOMP::PairLJCutCoulDSFOMP(LAMMPS *lmp) :
|
||||
@ -91,7 +86,7 @@ void PairLJCutCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double prefactor,erfcc,erfcd,t;
|
||||
double prefactor,erfcc,erfcd,grij;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -152,9 +147,9 @@ void PairLJCutCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (rsq < cut_coulsq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e*qtmp*q[j]/r;
|
||||
erfcd = exp(-alpha*alpha*r*r);
|
||||
t = 1.0 / (1.0 + EWALD_P*alpha*r);
|
||||
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
|
||||
grij = alpha*r;
|
||||
erfcd = expmsq(grij);
|
||||
erfcc = my_erfcx(grij) * erfcd;
|
||||
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
|
||||
r*f_shift) * r;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -18,18 +18,14 @@
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -91,7 +87,7 @@ void PairLJCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -143,18 +139,17 @@ void PairLJCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -149,7 +149,7 @@ void PairLJCutCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -187,7 +187,7 @@ void PairLJCutTholeLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qi*qj * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -26,14 +26,6 @@
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutTIP4PCutOMP::PairLJCutTIP4PCutOMP(LAMMPS *lmp) :
|
||||
|
||||
@ -20,19 +20,15 @@
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "error.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -140,7 +136,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
|
||||
double factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double v[6];
|
||||
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
|
||||
dbl3_t x1,x2,xH1,xH2;
|
||||
@ -304,11 +300,10 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (CTABLE || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
@ -317,7 +312,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -927,7 +927,7 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -18,12 +18,16 @@
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "lj_sdk_common.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
using namespace LJSDKParms;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -129,21 +133,12 @@ void PairLJSDKCoulLongOMP::eval_thr(int iifrom, int iito, ThrData * const thr)
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
const double A1 = 0.254829592;
|
||||
const double A2 = -0.284496736;
|
||||
const double A3 = 1.421413741;
|
||||
const double A4 = -1.453152027;
|
||||
const double A5 = 1.061405429;
|
||||
const double EWALD_F = 1.12837917;
|
||||
const double INV_EWALD_P = 1.0/0.3275911;
|
||||
|
||||
const double r = sqrt(rsq);
|
||||
const double grij = g_ewald * r;
|
||||
const double expm2 = exp(-grij*grij);
|
||||
const double t = INV_EWALD_P / (INV_EWALD_P + grij);
|
||||
const double erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
const double expm2 = expmsq(grij);
|
||||
const double erfc = my_erfcx(grij) * expm2;
|
||||
const double prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (EFLAG) ecoul = prefactor*erfc;
|
||||
if (sbindex) {
|
||||
const double adjust = (1.0-special_coul[sbindex])*prefactor;
|
||||
@ -154,7 +149,7 @@ void PairLJSDKCoulLongOMP::eval_thr(int iifrom, int iito, ThrData * const thr)
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
const double fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
const double table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
|
||||
|
||||
@ -152,7 +152,7 @@ void PairLJSDKCoulMSMOMP::eval_msm_thr(int iifrom, int iito, ThrData * const thr
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
const double fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
const double table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
|
||||
|
||||
@ -17,19 +17,15 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -91,7 +87,7 @@ void PairNMCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,factor_coul,factor_lj;
|
||||
double forcecoul,forcenm,rminv,rninv;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
int *ilist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -153,11 +149,10 @@ void PairNMCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (EFLAG) ecoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
@ -168,7 +163,7 @@ void PairNMCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG)
|
||||
|
||||
@ -166,7 +166,7 @@ void PairTableOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & tb->nmask;
|
||||
itable >>= tb->nshiftbits;
|
||||
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
|
||||
fraction = (rsq - tb->rsq[itable]) * tb->drsq[itable];
|
||||
value = tb->f[itable] + fraction*tb->df[itable];
|
||||
fpair = factor_lj * value;
|
||||
}
|
||||
|
||||
@ -26,14 +26,6 @@
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairTIP4PCutOMP::PairTIP4PCutOMP(LAMMPS *lmp) :
|
||||
|
||||
@ -21,18 +21,14 @@
|
||||
#include "neighbor.h"
|
||||
#include "error.h"
|
||||
#include "memory.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -140,7 +136,7 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,forcecoul,cforce;
|
||||
double factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double grij,expm2,prefactor,erfc;
|
||||
double v[6];
|
||||
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
|
||||
dbl3_t x1,x2,xH1,xH2;
|
||||
@ -278,11 +274,10 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
if (CTABLE || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
expm2 = expmsq(grij);
|
||||
erfc = my_erfcx(grij) * expm2;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
@ -291,7 +286,7 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
fraction = (rsq - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
|
||||
@ -27,20 +27,14 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulDSF::PairCoulDSF(LAMMPS *lmp) : Pair(lmp) {}
|
||||
@ -64,7 +58,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
|
||||
int i,j,ii,jj,inum,jnum;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double r,rsq,r2inv,forcecoul,factor_coul;
|
||||
double prefactor,erfcc,erfcd,t;
|
||||
double prefactor,erfcc,erfcd,grij;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
@ -115,9 +109,9 @@ void PairCoulDSF::compute(int eflag, int vflag)
|
||||
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e*qtmp*q[j]/r;
|
||||
erfcd = exp(-alpha*alpha*rsq);
|
||||
t = 1.0 / (1.0 + EWALD_P*alpha*r);
|
||||
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
|
||||
grij = alpha*r;
|
||||
erfcd = expmsq(grij);
|
||||
erfcc = my_erfcx(grij) * erfcd;
|
||||
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
|
||||
r*f_shift) * r;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
@ -295,7 +289,7 @@ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,erfcc,erfcd,prefactor,t;
|
||||
double r2inv,r,erfcc,erfcd,prefactor,grij;
|
||||
double forcecoul,phicoul;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
@ -303,15 +297,16 @@ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double eng = 0.0;
|
||||
if (rsq < cut_coulsq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = factor_coul * force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
erfcd = exp(-alpha*alpha*rsq);
|
||||
t = 1.0 / (1.0 + EWALD_P*alpha*r);
|
||||
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
|
||||
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
grij = alpha * r;
|
||||
erfcd = expmsq(grij);
|
||||
erfcc = my_erfcx(grij) * erfcd;
|
||||
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS*erfcd +
|
||||
r*f_shift) * r;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
|
||||
eng += phicoul;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
|
||||
@ -27,20 +27,14 @@
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "math_special.h"
|
||||
#include "math_const.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathSpecial;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulDSF::PairLJCutCoulDSF(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -77,7 +71,7 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag)
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double prefactor,erfcc,erfcd,t;
|
||||
double prefactor,erfcc,erfcd,grij;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -139,9 +133,9 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag)
|
||||
if (rsq < cut_coulsq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e*qtmp*q[j]/r;
|
||||
erfcd = exp(-alpha*alpha*r*r);
|
||||
t = 1.0 / (1.0 + EWALD_P*alpha*r);
|
||||
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
|
||||
grij = alpha*r;
|
||||
erfcd = expmsq(grij);
|
||||
erfcc = my_erfcx(grij) * erfcd;
|
||||
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
|
||||
r*f_shift) * r;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
@ -153,7 +153,7 @@ void PairTable::compute(int eflag, int vflag)
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & tb->nmask;
|
||||
itable >>= tb->nshiftbits;
|
||||
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
|
||||
fraction = (rsq - tb->rsq[itable]) * tb->drsq[itable];
|
||||
value = tb->f[itable] + fraction*tb->df[itable];
|
||||
fpair = factor_lj * value;
|
||||
}
|
||||
@ -1023,7 +1023,7 @@ double PairTable::single(int i, int j, int itype, int jtype, double rsq,
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & tb->nmask;
|
||||
itable >>= tb->nshiftbits;
|
||||
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
|
||||
fraction = (rsq - tb->rsq[itable]) * tb->drsq[itable];
|
||||
value = tb->f[itable] + fraction*tb->df[itable];
|
||||
fforce = factor_lj * value;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user