add tabulation for long-range coulomb

This commit is contained in:
Axel Kohlmeyer 2020-08-06 16:09:36 -04:00
parent 8826ea91e2
commit dae97e1151
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
1 changed files with 31 additions and 10 deletions

View File

@ -86,8 +86,9 @@ void PairBuckCoulLongOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
int i,j,ii,jj,jnum,itype,jtype,itable;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double fraction,table;
double rsq,r2inv,r6inv,r,rexp,forcecoul,forcebuck,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -139,6 +140,7 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
r = sqrt(rsq);
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -146,6 +148,20 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*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]) {
@ -167,7 +183,12 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
if (EFLAG) {
if (rsq < cut_coulsq) {
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]) {