add missing coulomb tabulation to pair style lj/class2/coul/long (9ed2824d) · Commits · 郑智淋 / lammps

src/USER-OMP/pair_lj_class2_coul_long_omp.cpp

+32 −10

Original line number	Diff line number	Diff line
		@@ -39,6 +39,7 @@ PairLJClass2CoulLongOMP::PairLJClass2CoulLongOMP(LAMMPS *lmp) :
		{
		suffix_flag \|= Suffix::OMP;
		respa_enable = 0;
		cut_respa = NULL;
		}

		/* ---------------------------------------------------------------------- */
		@@ -85,8 +86,9 @@ void PairLJClass2CoulLongOMP::compute(int eflag, int vflag)
		template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
		void PairLJClass2CoulLongOMP::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 r,rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
		double grij,expm2,prefactor,t,erfc;
		int ilist,jlist,numneigh,*firstneigh;
		@@ -137,6 +139,7 @@ void PairLJClass2CoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
		r2inv = 1.0/rsq;

		if (rsq < cut_coulsq) {
		if (!ncoultablebits \|\| rsq <= tabinnersq) {
		r = sqrt(rsq);
		grij = g_ewald * r;
		expm2 = exp(-grij*grij);
		@@ -145,6 +148,20 @@ void PairLJClass2CoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
		prefactor = qqrd2e * qtmp*q[j]/r;
		forcecoul = prefactor * (erfc + EWALD_Fgrijexpm2);
		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 = qtmpq[j] table;
		if (factor_coul < 1.0) {
		table = ctable[itable] + fraction*dctable[itable];
		prefactor = qtmpq[j] table;
		forcecoul -= (1.0-factor_coul)*prefactor;
		}
		}
		} else forcecoul = 0.0;

		if (rsq < cut_ljsq[itype][jtype]) {
		@@ -168,7 +185,12 @@ void PairLJClass2CoulLongOMP::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 = qtmpq[j] table;
		}
		if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
		} else ecoul = 0.0;

Admin message