Compare commits
25 Commits
patch_28Ma
...
coulomb-re
| Author | SHA1 | Date | |
|---|---|---|---|
| a31134a16f | |||
| ae56b9ad89 | |||
| 4466d9fb4a | |||
| ac1aa9edea | |||
| c733204a70 | |||
| a89246fbfc | |||
| f739a94eeb | |||
| 9d8d88d020 | |||
| 6410797697 | |||
| 0898ffe336 | |||
| d84ffc5309 | |||
| 124bd99751 | |||
| 85a24bda21 | |||
| d759949d03 | |||
| 25904e80f0 | |||
| cfb2353313 | |||
| bd49f4d665 | |||
| b2ef952187 | |||
| 381d87b79b | |||
| 8778452eaf | |||
| 6ae15d0801 | |||
| 2826449b52 | |||
| c7c6c44add | |||
| 8e145e2b9f | |||
| 0f76bbc718 |
@ -1,7 +1,7 @@
|
||||
<!-- HTML_ONLY -->
|
||||
<HEAD>
|
||||
<TITLE>LAMMPS Users Manual</TITLE>
|
||||
<META NAME="docnumber" CONTENT="28 Mar 2017 version">
|
||||
<META NAME="docnumber" CONTENT="31 Mar 2017 version">
|
||||
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
|
||||
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
|
||||
</HEAD>
|
||||
@ -21,7 +21,7 @@
|
||||
<H1></H1>
|
||||
|
||||
LAMMPS Documentation :c,h3
|
||||
28 Mar 2017 version :c,h4
|
||||
31 Mar 2017 version :c,h4
|
||||
|
||||
Version info: :h4
|
||||
|
||||
|
||||
@ -1686,7 +1686,7 @@ nph) and Berendsen:
|
||||
The "fix npt"_fix_nh.html commands include a Nose-Hoover thermostat
|
||||
and barostat. "Fix nph"_fix_nh.html is just a Nose/Hoover barostat;
|
||||
it does no thermostatting. Both "fix nph"_fix_nh.html and "fix
|
||||
press/bernendsen"_fix_press_berendsen.html can be used in conjunction
|
||||
press/berendsen"_fix_press_berendsen.html can be used in conjunction
|
||||
with any of the thermostatting fixes.
|
||||
|
||||
As with the thermostats, "fix npt"_fix_nh.html and "fix
|
||||
|
||||
@ -90,8 +90,8 @@ the difference between the {charmm} and {charmmfsh} styles is in the
|
||||
computation of the 1-4 non-bond interactions, though only if the
|
||||
distance between the two atoms is within the switching region of the
|
||||
pairwise potential defined by the corresponding CHARMM pair style,
|
||||
i.e. between the inner and outer cutoffs specified for the pair style.
|
||||
The {charmmfsh} style should only be used when using the "pair_style
|
||||
i.e. within the outer cutoff specified for the pair style. The
|
||||
{charmmfsh} style should only be used when using the "pair_style
|
||||
lj/charmmfsw/coul/charmmfsh"_pair_charmm.html to make the Coulombic
|
||||
pairwise calculations consistent. Use the {charmm} style with
|
||||
long-range Coulombics or the older "pair_style
|
||||
|
||||
@ -464,6 +464,7 @@ pair_nb3b_harmonic.html
|
||||
pair_nm.html
|
||||
pair_none.html
|
||||
pair_oxdna.html
|
||||
pair_oxdna2.html
|
||||
pair_peri.html
|
||||
pair_polymorphic.html
|
||||
pair_quip.html
|
||||
|
||||
@ -49,8 +49,8 @@ args = list of arguments for a particular style :ul
|
||||
|
||||
pair_style lj/charmm/coul/charmm 8.0 10.0
|
||||
pair_style lj/charmm/coul/charmm 8.0 10.0 7.0 9.0
|
||||
pair_style lj/charmmfsw/coul/charmmfsh 8.0 10.0
|
||||
pair_style lj/charmmfsw/coul/charmmfsh 8.0 10.0 7.0 9.0
|
||||
pair_style lj/charmmfsw/coul/charmmfsh 10.0 12.0
|
||||
pair_style lj/charmmfsw/coul/charmmfsh 10.0 12.0 9.0
|
||||
pair_coeff * * 100.0 2.0
|
||||
pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
|
||||
|
||||
|
||||
@ -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];
|
||||
|
||||
@ -29,19 +29,29 @@ action () {
|
||||
# list of files with dependcies
|
||||
|
||||
action bond_oxdna_fene.cpp bond_fene.h
|
||||
action bond_oxdna2_fene.cpp bond_fene.h
|
||||
action bond_oxdna_fene.h bond_fene.h
|
||||
action bond_oxdna2_fene.h bond_fene.h
|
||||
action fix_nve_dotc_langevin.cpp atom_vec_ellipsoid.h
|
||||
action fix_nve_dotc_langevin.h atom_vec_ellipsoid.h
|
||||
action fix_nve_dot.cpp atom_vec_ellipsoid.h
|
||||
action fix_nve_dot.h atom_vec_ellipsoid.h
|
||||
action mf_oxdna.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna_coaxstk.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_coaxstk.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna_coaxstk.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_coaxstk.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna_excv.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_excv.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna_excv.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_excv.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna_hbond.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna_hbond.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna_stk.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_stk.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna_stk.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_stk.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna_xstk.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna_xstk.h atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_dh.cpp atom_vec_ellipsoid.h
|
||||
action pair_oxdna2_dh.h atom_vec_ellipsoid.h
|
||||
|
||||
@ -2,9 +2,9 @@ This package contains a LAMMPS implementation of coarse-grained
|
||||
models of DNA, which can be used to model sequence-specific
|
||||
DNA strands.
|
||||
|
||||
See the doc pages and [1,2] for the individual bond and pair styles.
|
||||
See the doc pages and [1,2,3] for the individual bond and pair styles.
|
||||
The packages contains also a new Langevin-type rigid-body integrator,
|
||||
which has also its own doc page and is explained in [3].
|
||||
which has also its own doc page and is explained in [4].
|
||||
|
||||
[1] T. Ouldridge, A. Louis, J. Doye, "Structural, mechanical,
|
||||
and thermodynamic properties of a coarse-grained DNA model",
|
||||
@ -13,16 +13,20 @@ J. Chem. Phys. 134, 085101 (2011).
|
||||
[2] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA
|
||||
self-assembly, DPhil. University of Oxford (2011).
|
||||
|
||||
[3] R. Davidchack, T. Ouldridge, M. Tretyakov, "New Langevin and
|
||||
[3] B.E. Snodin, F. Randisi, M. Mosayebi, et al., Introducing
|
||||
Improved Structural Properties and Salt Dependence into a Coarse-Grained
|
||||
Model of DNA, J. Chem. Phys. 142, 234901 (2015).
|
||||
|
||||
[4] R. Davidchack, T. Ouldridge, M. Tretyakov, "New Langevin and
|
||||
gradient thermostats for rigid body dynamics", J. Chem. Phys. 142,
|
||||
144114 (2015).
|
||||
|
||||
Example input and data files can be found in
|
||||
/examples/USER/cgdna/examples/duplex1/ and /duplex2/. A simple python
|
||||
setup tool which creates single straight or helical DNA strands as
|
||||
well as DNA duplexes and arrays of duplexes can be found in
|
||||
/examples/USER/cgdna/util/. A technical report with more information
|
||||
on the model, the structure of the input and data file, the setup tool
|
||||
/examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/. Python setup
|
||||
tools which create single straight or helical DNA strands as
|
||||
well as DNA duplexes or arrays of duplexes can be found in
|
||||
/examples/USER/cgdna/util/. A technical report with more information
|
||||
on the models, the structure of the input and data file, the setup tool
|
||||
and the performance of the LAMMPS-implementation of oxDNA can be found
|
||||
in /doc/src/PDF/USER-CGDNA-overview.pdf.
|
||||
|
||||
@ -35,6 +39,7 @@ of the LAMMPS manual).
|
||||
The creator of this package is:
|
||||
|
||||
Dr Oliver Henrich
|
||||
University of Strathclyde, UK
|
||||
University of Edinburgh, UK
|
||||
ohenrich@ph.ed.ac.uk
|
||||
o.henrich@epcc.ed.ac.uk
|
||||
@ -45,6 +50,8 @@ o.henrich@epcc.ed.ac.uk
|
||||
|
||||
bond_oxdna_fene.cpp: backbone connectivity, a modified FENE potential
|
||||
|
||||
bond_oxdna2_fene.cpp: corresponding bond style in oxDNA2 (see [3])
|
||||
|
||||
** Pair styles provided by this package:
|
||||
|
||||
pair_oxdna_excv.cpp: excluded volume interaction between the nucleotides
|
||||
@ -59,6 +66,13 @@ pair_oxdna_xstk.cpp: cross-stacking interaction between nucleotides
|
||||
|
||||
pair_oxdna_coaxstk.cpp: coaxial stacking interaction between nucleotides
|
||||
|
||||
|
||||
pair_oxdna2_excv.cpp, pair_oxdna2_stk.cpp, pair_oxdna2_coaxstk.cpp:
|
||||
corresponding pair styles in oxDNA2 (see [3])
|
||||
|
||||
pair_oxdna2_dh.cpp: Debye-Hueckel electrostatic interaction between backbone
|
||||
sites
|
||||
|
||||
** Fixes provided by this package:
|
||||
|
||||
fix_nve_dotc_langevin.cpp: fix for Langevin-type rigid body integrator "C"
|
||||
|
||||
@ -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) {
|
||||
|
||||
@ -254,10 +254,9 @@ double FixHalt::tlimit()
|
||||
bigint final = update->firststep +
|
||||
static_cast<bigint> (tratio*value/cpu * elapsed);
|
||||
nextstep = (final/nevery)*nevery + nevery;
|
||||
if (nextstep == update->ntimestep) nextstep += nevery;
|
||||
tratio = 1.0;
|
||||
}
|
||||
|
||||
//printf("EVAL %ld %g %d\n",update->ntimestep,cpu,nevery);
|
||||
|
||||
return cpu;
|
||||
}
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
@ -1 +1 @@
|
||||
#define LAMMPS_VERSION "28 Mar 2017"
|
||||
#define LAMMPS_VERSION "31 Mar 2017"
|
||||
|
||||
Reference in New Issue
Block a user