diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp index 0d0a9be416..ba24fe1864 100644 --- a/src/CLASS2/pair_lj_class2_coul_long.cpp +++ b/src/CLASS2/pair_lj_class2_coul_long.cpp @@ -42,6 +42,7 @@ using namespace MathConst; PairLJClass2CoulLong::PairLJClass2CoulLong(LAMMPS *lmp) : Pair(lmp) { ewaldflag = pppmflag = 1; + ftable = NULL; } /* ---------------------------------------------------------------------- */ @@ -68,8 +69,9 @@ PairLJClass2CoulLong::~PairLJClass2CoulLong() void PairLJClass2CoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itable,itype,jtype; 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 factor_coul,factor_lj; @@ -122,14 +124,29 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) r2inv = 1.0/rsq; 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; - 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 = 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; + } 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]) { @@ -152,7 +169,12 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) 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]) { @@ -270,9 +292,12 @@ void PairLJClass2CoulLong::init_style() // insure use of KSpace long-range solver, set g_ewald - if (force->kspace == NULL) + if (force->kspace == NULL) error->all(FLERR,"Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; + + // setup force tables + if (ncoultablebits) init_tables(cut_coul,NULL); } /* ---------------------------------------------------------------------- @@ -400,6 +425,8 @@ void PairLJClass2CoulLong::write_restart_settings(FILE *fp) fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&ncoultablebits,sizeof(int),1,fp); + fwrite(&tabinner,sizeof(double),1,fp); } /* ---------------------------------------------------------------------- @@ -414,12 +441,16 @@ void PairLJClass2CoulLong::read_restart_settings(FILE *fp) fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); fread(&tail_flag,sizeof(int),1,fp); + fread(&ncoultablebits,sizeof(int),1,fp); + fread(&tabinner,sizeof(double),1,fp); } MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); + MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ @@ -430,18 +461,34 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, double &fforce) { double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor; - double forcecoul,forcelj,phicoul,philj; + double fraction,table,forcecoul,forcelj,phicoul,philj; + int itable; r2inv = 1.0/rsq; 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; - prefactor = force->qqrd2e * atom->q[i]*atom->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 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->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 = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } } else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); @@ -453,7 +500,12 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, double eng = 0.0; if (rsq < cut_coulsq) { - phicoul = prefactor*erfc; + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; eng += phicoul; }