git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8159 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp

@@ -175,12 +175,10 @@ void PairLJSDKCoulLong::eval()
 	    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 (EFLAG)
+	    if (EFLAG) ecoul = prefactor*erfc;
-	      ecoul = prefactor*erfc;
 	    if (factor_coul < 1.0) {
 	      forcecoul -= (1.0-factor_coul)*prefactor;
-	      if (EFLAG)
+	      if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
-		ecoul -= (1.0-factor_coul)*prefactor;
 	    }
 	  } else {
 	    union_int_float_t rsq_lookup;
@@ -190,14 +188,12 @@ void PairLJSDKCoulLong::eval()
 	    fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
 	    table = ftable[itable] + fraction*dftable[itable];
 	    forcecoul = qtmp*q[j] * table;
-	    if (EFLAG)
+	    if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
-	      ecoul = 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;
-	      if (EFLAG)
+	      if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
-		ecoul -= (1.0-factor_coul)*prefactor;
 	    }
 	  }
 	}
@@ -231,8 +227,7 @@ void PairLJSDKCoulLong::eval()
 			     - lj4[itype][jtype]) - offset[itype][jtype];
 	  }
 	  forcelj *= factor_lj;
-	  if (EFLAG)
+	  if (EFLAG) evdwl *= factor_lj;
-	    evdwl *= factor_lj;
 	}
 
 	fpair = (forcecoul + forcelj) * r2inv;
@@ -248,7 +243,6 @@ void PairLJSDKCoulLong::eval()
 
 	if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
 			     evdwl,ecoul,fpair,delx,dely,delz);
-
       }
     }
     f[i][0] += fxtmp;

src/USER-OMP/improper_class2_omp.cpp

@@ -79,7 +79,7 @@ void ImproperClass2OMP::compute(int eflag, int vflag)
   } // end of omp parallel region
 }
 
-template <const int EVFLAG, const int EFLAG, const int NEWTON_BOND>
+template <int EVFLAG, int EFLAG, int NEWTON_BOND>
 void ImproperClass2OMP::eval(int nfrom, int nto, ThrData * const thr)
 {
   int i1,i2,i3,i4,i,j,k,n,type;

src/USER-OMP/pair_buck_coul_omp.cpp

@@ -39,6 +39,7 @@ PairBuckCoulOMP::PairBuckCoulOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_coul_long_omp.cpp

@@ -38,6 +38,7 @@ PairCoulLongOMP::PairCoulLongOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_lj96_cut_omp.cpp

@@ -30,6 +30,7 @@ PairLJ96CutOMP::PairLJ96CutOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_lj_charmm_coul_long_omp.cpp

@@ -23,14 +23,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
-
 /* ---------------------------------------------------------------------- */
 
 PairLJCharmmCoulLongOMP::PairLJCharmmCoulLongOMP(LAMMPS *lmp) :
@@ -38,6 +30,7 @@ PairLJCharmmCoulLongOMP::PairLJCharmmCoulLongOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -84,101 +77,121 @@ void PairLJCharmmCoulLongOMP::compute(int eflag, int vflag)
 template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
 void PairLJCharmmCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 {
-  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,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
-  double grij,expm2,prefactor,t,erfc;
-  double philj,switch1,switch2;
-  int *ilist,*jlist,*numneigh,**firstneigh;
-
-  evdwl = ecoul = 0.0;
 
   const double * const * const x = atom->x;
   double * const * const f = thr->get_f();
   const double * const q = atom->q;
   const int * const type = atom->type;
-  const int nlocal = atom->nlocal;
   const double * const special_coul = force->special_coul;
   const double * const special_lj = force->special_lj;
   const double qqrd2e = force->qqrd2e;
-  double fxtmp,fytmp,fztmp;
+  const double inv_denom_lj = 1.0/denom_lj;
 
-  ilist = list->ilist;
+  const int * const ilist = list->ilist;
-  numneigh = list->numneigh;
+  const int * const numneigh = list->numneigh;
-  firstneigh = list->firstneigh;
+  const int * const * const firstneigh = list->firstneigh;
+  const int nlocal = atom->nlocal;
 
   // loop over neighbors of my atoms
 
-  for (ii = iifrom; ii < iito; ++ii) {
+  for (int ii = iifrom; ii < iito; ++ii) {
 
-    i = ilist[ii];
+    const int i = ilist[ii];
-    qtmp = q[i];
+    const int itype = type[i];
-    xtmp = x[i][0];
+    const double qtmp = q[i];
-    ytmp = x[i][1];
+    const double xtmp = x[i][0];
-    ztmp = x[i][2];
+    const double ytmp = x[i][1];
-    itype = type[i];
+    const double ztmp = x[i][2];
-    jlist = firstneigh[i];
+    double fxtmp,fytmp,fztmp;
-    jnum = numneigh[i];
     fxtmp=fytmp=fztmp=0.0;
 
-    for (jj = 0; jj < jnum; jj++) {
+    const int * const jlist = firstneigh[i];
-      j = jlist[jj];
+    const int jnum = numneigh[i];
-      factor_lj = special_lj[sbmask(j)];
-      factor_coul = special_coul[sbmask(j)];
-      j &= NEIGHMASK;
 
-      delx = xtmp - x[j][0];
+    for (int jj = 0; jj < jnum; jj++) {
-      dely = ytmp - x[j][1];
+      double forcecoul, forcelj, evdwl, ecoul;
-      delz = ztmp - x[j][2];
+      forcecoul = forcelj = evdwl = ecoul = 0.0;
-      rsq = delx*delx + dely*dely + delz*delz;
+
-      jtype = type[j];
+      const int sbindex = sbmask(jlist[jj]);
+      const int j = jlist[jj] & NEIGHMASK;
+
+      const double delx = xtmp - x[j][0];
+      const double dely = ytmp - x[j][1];
+      const double delz = ztmp - x[j][2];
+      const double rsq = delx*delx + dely*dely + delz*delz;
+      const int jtype = type[j];
 
       if (rsq < cutsq[itype][jtype]) {
-	r2inv = 1.0/rsq;
+	const double r2inv = 1.0/rsq;
 
 	if (rsq < cut_coulsq) {
 	  if (!ncoultablebits || rsq <= tabinnersq) {
-	    r = sqrt(rsq);
+	    const double A1 =  0.254829592;
-	    grij = g_ewald * r;
+	    const double A2 = -0.284496736;
-	    expm2 = exp(-grij*grij);
+	    const double A3 =  1.421413741;
-	    t = 1.0 / (1.0 + EWALD_P*grij);
+	    const double A4 = -1.453152027;
-	    erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+	    const double A5 =  1.061405429;
-	    prefactor = qqrd2e * qtmp*q[j]/r;
+	    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 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 (EFLAG) ecoul = prefactor*erfc;
+	    if (sbindex) {
+	      const double adjust = (1.0-special_coul[sbindex])*prefactor;
+	      forcecoul -= adjust;
+	      if (EFLAG) ecoul -= adjust;
+	    }
 	  } else {
 	    union_int_float_t rsq_lookup;
 	    rsq_lookup.f = rsq;
-	    itable = rsq_lookup.i & ncoulmask;
+	    const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
-	    itable >>= ncoulshiftbits;
+	    const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
-	    fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+	    const double table = ftable[itable] + fraction*dftable[itable];
-	    table = ftable[itable] + fraction*dftable[itable];
 	    forcecoul = qtmp*q[j] * table;
-	    if (factor_coul < 1.0) {
+	    if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
-	      table = ctable[itable] + fraction*dctable[itable];
+	    if (sbindex) {
-	      prefactor = qtmp*q[j] * table;
+	      const double table2 = ctable[itable] + fraction*dctable[itable];
-	      forcecoul -= (1.0-factor_coul)*prefactor;
+	      const double prefactor = qtmp*q[j] * table2;
+	      const double adjust = (1.0-special_coul[sbindex])*prefactor;
+	      forcecoul -= adjust;
+	      if (EFLAG) ecoul -= adjust;
+	    }
 	  }
 	}
-	} else forcecoul = 0.0;
 
 	if (rsq < cut_ljsq) {
-	  r6inv = r2inv*r2inv*r2inv;
+	  const double r6inv = r2inv*r2inv*r2inv;
-	  jtype = type[j];
 	  forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+	  if (EFLAG) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+
 	  if (rsq > cut_lj_innersq) {
-	    switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+	    const double drsq = cut_ljsq - rsq;
-	      (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+	    const double cut2 = (rsq - cut_lj_innersq) * drsq;
-	    switch2 = 12.0*rsq * (cut_ljsq-rsq) * 
+	    const double switch1 = drsq * (drsq*drsq + 3.0*cut2) * inv_denom_lj;
-	      (rsq-cut_lj_innersq) / denom_lj;
+	    const double switch2 = 12.0*rsq * cut2 * inv_denom_lj;
-	    philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
+	    if (EFLAG) {
+	      forcelj = forcelj*switch1 + evdwl*switch2;
+	      evdwl *= switch1;
+	    } else {
+	      const double philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
 	      forcelj =  forcelj*switch1 + philj*switch2;
 	    }
-	  forcelj *= factor_lj;
+	  }
-	} else forcelj = 0.0;
 
-	fpair = (forcecoul + forcelj) * r2inv;
+	  if (sbindex) {
+	    const double factor_lj = special_lj[sbindex];
+	    forcelj *= factor_lj;
+	    if (EFLAG) evdwl *= factor_lj;
+	  }
+
+	}
+	const double fpair = (forcecoul + forcelj) * r2inv;
 
 	fxtmp += delx*fpair;
 	fytmp += dely*fpair;
@@ -189,29 +202,7 @@ void PairLJCharmmCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 	  f[j][2] -= delz*fpair;
 	}
 
-	if (EFLAG) {
+	if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
-	  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) {
-	    evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
-	    if (rsq > cut_lj_innersq) {
-	      switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
-		(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
-	      evdwl *= switch1;
-	    }
-	    evdwl *= factor_lj;
-	  } else evdwl = 0.0;
-	}
-	
-	if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
 				 evdwl,ecoul,fpair,delx,dely,delz,thr);
       }
     }

src/USER-OMP/pair_lj_coul_omp.cpp

@@ -39,6 +39,7 @@ PairLJCoulOMP::PairLJCoulOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_lj_cut_coul_long_omp.cpp

@@ -38,6 +38,7 @@ PairLJCutCoulLongOMP::PairLJCutCoulLongOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp

@@ -145,7 +145,7 @@ void PairLJCutCoulLongTIP4POMP::compute(int eflag, int vflag)
 
 /* ---------------------------------------------------------------------- */
 
-template <const int CTABLE, const int EVFLAG, const int EFLAG, const int VFLAG>
+template <int CTABLE, int EVFLAG, int EFLAG, int VFLAG>
 void PairLJCutCoulLongTIP4POMP::eval(int iifrom, int iito, ThrData * const thr)
 {
   int i,j,ii,jj,jnum,itype,jtype,itable;

src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp

@@ -167,7 +167,7 @@ void PairLJCutCoulPPPMTIP4POMP::compute(int eflag, int vflag)
 
 /* ---------------------------------------------------------------------- */
 
-template <const int CTABLE, const int EVFLAG, const int EFLAG, const int VFLAG>
+template <int CTABLE, int EVFLAG, int EFLAG, int VFLAG>
 void PairLJCutCoulPPPMTIP4POMP::eval(int iifrom, int iito, ThrData * const thr)
 {
   int i,j,ii,jj,jnum,itype,jtype,itable;

src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.h

@@ -48,7 +48,7 @@ class PairLJCutCoulPPPMTIP4POMP : public PairLJCutCoulLongTIP4P, public ThrOMP {
   void find_M_permissive(int, int &, int &, double *);
 
  private:
-  template <const int, const int, const int, const int>
+  template <int, int, int, int>
   void eval(int ifrom, int ito, ThrData * const thr);
 
   class PPPMTIP4PProxy *kspace;

src/USER-OMP/pair_lj_cut_omp.cpp

@@ -30,6 +30,7 @@ PairLJCutOMP::PairLJCutOMP(LAMMPS *lmp) :
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;
+  cut_respa = NULL;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_lj_sdk_coul_long_omp.cpp

@@ -25,15 +25,6 @@
 #include "suffix.h"
 using namespace LAMMPS_NS;
 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
-
 /* ---------------------------------------------------------------------- */
 
 PairLJSDKCoulLongOMP::PairLJSDKCoulLongOMP(LAMMPS *lmp) :
@@ -87,91 +78,90 @@ void PairLJSDKCoulLongOMP::compute(int eflag, int vflag)
 template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
 void PairLJSDKCoulLongOMP::eval_thr(int iifrom, int iito, ThrData * const thr)
 {
-  int i,ii,j,jj,jtype,itable;
-  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;
 
   const double * const * const x = atom->x;
   double * const * const f = thr->get_f();
   const double * const q = atom->q;
   const int * const type = atom->type;
-  const int nlocal = atom->nlocal;
   const double * const special_coul = force->special_coul;
   const double * const special_lj = force->special_lj;
   const double qqrd2e = force->qqrd2e;
-  double fxtmp,fytmp,fztmp;
 
   const int * const ilist = list->ilist;
   const int * const numneigh = list->numneigh;
   const int * const * const firstneigh = list->firstneigh;
+  const int nlocal = atom->nlocal;
 
   // loop over neighbors of my atoms
 
-  for (ii = iifrom; ii < iito; ++ii) {
+  for (int ii = iifrom; ii < iito; ++ii) {
 
-    i = ilist[ii];
+    const int i = ilist[ii];
-    qtmp = q[i];
+    const int itype = type[i];
-    xtmp = x[i][0];
+    const double qtmp = q[i];
-    ytmp = x[i][1];
+    const double xtmp = x[i][0];
-    ztmp = x[i][2];
+    const double ytmp = x[i][1];
+    const double ztmp = x[i][2];
+    double fxtmp,fytmp,fztmp;
     fxtmp=fytmp=fztmp=0.0;
 
-    const int itype = type[i];
     const int * const jlist = firstneigh[i];
     const int jnum = numneigh[i];
 
-    for (jj = 0; jj < jnum; jj++) {
+    for (int jj = 0; jj < jnum; jj++) {
+      double forcecoul, forcelj, evdwl, ecoul;
       forcecoul = forcelj = evdwl = ecoul = 0.0;
 
-      j = jlist[jj];
+      const int sbindex = sbmask(jlist[jj]);
-      factor_lj = special_lj[sbmask(j)];
+      const int j = jlist[jj] & NEIGHMASK;
-      factor_coul = special_coul[sbmask(j)];
-      j &= NEIGHMASK;
 
-      delx = xtmp - x[j][0];
+      const double delx = xtmp - x[j][0];
-      dely = ytmp - x[j][1];
+      const double dely = ytmp - x[j][1];
-      delz = ztmp - x[j][2];
+      const double delz = ztmp - x[j][2];
-      rsq = delx*delx + dely*dely + delz*delz;
+      const double rsq = delx*delx + dely*dely + delz*delz;
-      jtype = type[j];
+      const int jtype = type[j];
 
       if (rsq < cutsq[itype][jtype]) {
-	r2inv = 1.0/rsq;
+	const double r2inv = 1.0/rsq;
 	const int ljt = lj_type[itype][jtype];
 
 	if (rsq < cut_coulsq) {
 	  if (!ncoultablebits || rsq <= tabinnersq) {
-	    r = sqrt(rsq);
+	    const double A1 =  0.254829592;
-	    grij = g_ewald * r;
+	    const double A2 = -0.284496736;
-	    expm2 = exp(-grij*grij);
+	    const double A3 =  1.421413741;
-	    t = 1.0 / (1.0 + EWALD_P*grij);
+	    const double A4 = -1.453152027;
-	    erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+	    const double A5 =  1.061405429;
-	    prefactor = qqrd2e * qtmp*q[j]/r;
+	    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 prefactor = qqrd2e * qtmp*q[j]/r;
 	    forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
-	    if (EFLAG)
+	    if (EFLAG) ecoul = prefactor*erfc;
-	      ecoul = prefactor*erfc;
+	    if (sbindex) {
-	    if (factor_coul < 1.0) {
+	      const double adjust = (1.0-special_coul[sbindex])*prefactor;
-	      forcecoul -= (1.0-factor_coul)*prefactor;
+	      forcecoul -= adjust;
-	      if (EFLAG)
+	      if (EFLAG) ecoul -= adjust;
-		ecoul -= (1.0-factor_coul)*prefactor;
 	    }
 	  } else {
 	    union_int_float_t rsq_lookup;
 	    rsq_lookup.f = rsq;
-	    itable = rsq_lookup.i & ncoulmask;
+	    const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
-	    itable >>= ncoulshiftbits;
+	    const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
-	    fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+	    const double table = ftable[itable] + fraction*dftable[itable];
-	    table = ftable[itable] + fraction*dftable[itable];
 	    forcecoul = qtmp*q[j] * table;
-	    if (EFLAG)
+	    if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
-	      ecoul = qtmp*q[j] * table;
+	    if (sbindex) {
-	    if (factor_coul < 1.0) {
+	      const double table2 = ctable[itable] + fraction*dctable[itable];
-	      table = ctable[itable] + fraction*dctable[itable];
+	      const double prefactor = qtmp*q[j] * table2;
-	      prefactor = qtmp*q[j] * table;
+	      const double adjust = (1.0-special_coul[sbindex])*prefactor;
-	      forcecoul -= (1.0-factor_coul)*prefactor;
+	      forcecoul -= adjust;
-	      if (EFLAG)
+	      if (EFLAG) ecoul -= adjust;
-		ecoul -= (1.0-factor_coul)*prefactor;
 	    }
 	  }
 	}
@@ -204,12 +194,16 @@ void PairLJSDKCoulLongOMP::eval_thr(int iifrom, int iito, ThrData * const thr)
 	      evdwl = r6inv*(lj3[itype][jtype]*r6inv
 			     - lj4[itype][jtype]) - offset[itype][jtype];
 	  }
+
+	  if (sbindex) {
+	    const double factor_lj = special_lj[sbindex];
 	    forcelj *= factor_lj;
-	  if (EFLAG)
+	    if (EFLAG) evdwl *= factor_lj;
-	    evdwl *= factor_lj;
 	  }
 
-	fpair = (forcecoul + forcelj) * r2inv;
+	}
+
+	const double fpair = (forcecoul + forcelj) * r2inv;
 
 	fxtmp += delx*fpair;
 	fytmp += dely*fpair;
@@ -222,7 +216,6 @@ void PairLJSDKCoulLongOMP::eval_thr(int iifrom, int iito, ThrData * const thr)
 
 	if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
 				 evdwl,ecoul,fpair,delx,dely,delz,thr);
-
       }
     }
     f[i][0] += fxtmp;

src/XTC/dump_xtc.cpp

@@ -189,21 +189,6 @@ void DumpXTC::write_header(bigint nbig)
 
 /* ---------------------------------------------------------------------- */
 
-int DumpXTC::count()
-{
-  if (igroup == 0) return atom->nlocal;
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  int m = 0;
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) m++;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void DumpXTC::pack(int *ids)
 {
   int m,n;

src/XTC/dump_xtc.h

@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -49,7 +49,6 @@ class DumpXTC : public Dump {
   int modify_param(int, char **);
   void openfile();
   void write_header(bigint);
-  int count();
   void pack(int *);
   void write_data(int, double *);
   bigint memory_usage();

src/dump.cpp

@@ -66,6 +66,7 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   flush_flag = 1;
   format = NULL;
   format_user = NULL;
+  format_default = NULL;
   clearstep = 0;
   sort_flag = 0;
   append_flag = 0;
@@ -223,6 +224,21 @@ void Dump::init()
 
 /* ---------------------------------------------------------------------- */
 
+int Dump::count()
+{
+  if (igroup == 0) return atom->nlocal;
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int m = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) m++;
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Dump::write()
 {
   // if file per timestep, open new file

src/dump.h

@@ -102,7 +102,7 @@ class Dump : protected Pointers {
   virtual void openfile();
   virtual int modify_param(int, char **) {return 0;}
   virtual void write_header(bigint) = 0;
-  virtual int count() = 0;
+  virtual int count();
   virtual void pack(int *) = 0;
   virtual void write_data(int, double *) = 0;
 

src/dump_atom.cpp

@@ -145,21 +145,6 @@ void DumpAtom::write_header(bigint ndump)
 
 /* ---------------------------------------------------------------------- */
 
-int DumpAtom::count()
-{
-  if (igroup == 0) return atom->nlocal;
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  int m = 0;
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) m++;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void DumpAtom::pack(int *ids)
 {
   (this->*pack_choice)(ids);

src/dump_atom.h

@@ -37,7 +37,6 @@ class DumpAtom : public Dump {
   void init_style();
   int modify_param(int, char **);
   void write_header(bigint);
-  int count();
   void pack(int *);
   void write_data(int, double *);
 

src/dump_dcd.cpp

@@ -174,21 +174,6 @@ void DumpDCD::write_header(bigint n)
 
 /* ---------------------------------------------------------------------- */
 
-int DumpDCD::count()
-{
-  if (igroup == 0) return atom->nlocal;
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  int m = 0;
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) m++;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void DumpDCD::pack(int *ids)
 {
   int m,n;

src/dump_dcd.h

@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -20,9 +20,7 @@ DumpStyle(dcd,DumpDCD)
 #ifndef LMP_DUMP_DCD_H
 #define LMP_DUMP_DCD_H
 
-#include "stdio.h"
 #include "dump.h"
-#include "inttypes.h"
 
 namespace LAMMPS_NS {
 
@@ -41,7 +39,6 @@ class DumpDCD : public Dump {
   void init_style();
   void openfile();
   void write_header(bigint);
-  int count();
   void pack(int *);
   void write_data(int, double *);
   int modify_param(int, char **);

src/dump_xyz.cpp

@@ -17,6 +17,7 @@
 #include "group.h"
 #include "error.h"
 #include "memory.h"
+#include "update.h"
 
 using namespace LAMMPS_NS;
 
@@ -31,10 +32,30 @@ DumpXYZ::DumpXYZ(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
   sort_flag = 1;
   sortcol = 0;
 
-  char *str = (char *) "%d %g %g %g";
+  if (format_default) delete [] format_default;
+
+  char *str = (char *) "%s %g %g %g";
   int n = strlen(str) + 1;
   format_default = new char[n];
   strcpy(format_default,str);
+
+  ntypes = atom->ntypes;
+  typenames = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+DumpXYZ::~DumpXYZ()
+{
+  delete[] format_default;
+  format_default = NULL;
+
+  if (typenames) {
+    for (int i = 1; i <= ntypes; i++)
+      delete [] typenames[i];
+    delete [] typenames;
+    typenames = NULL;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -51,6 +72,17 @@ void DumpXYZ::init_style()
   strcpy(format,str);
   strcat(format,"\n");
 
+  // initialize typenames array to be backward compatible by default
+  // a 32-bit int can be maximally 10 digits plus sign
+
+  if (typenames == NULL) {
+    typenames = new char*[ntypes+1];
+    for (int itype = 1; itype <= ntypes; itype++) {
+      typenames[itype] = new char[12];
+      sprintf(typenames[itype],"%d",itype);
+    }
+  }
+
   // open single file, one time only
 
   if (multifile == 0) openfile();
@@ -58,31 +90,45 @@ void DumpXYZ::init_style()
 
 /* ---------------------------------------------------------------------- */
 
+int DumpXYZ::modify_param(int narg, char **arg)
+{
+  if (strcmp(arg[0],"element") == 0) {
+    if (narg < ntypes+1)
+      error->all(FLERR, "Dump modify element names do not match atom types");
+
+    if (typenames) {
+      for (int i = 1; i <= ntypes; i++)
+	delete [] typenames[i];
+
+      delete [] typenames;
+      typenames = NULL;
+    }
+
+    typenames = new char*[ntypes+1];
+    for (int itype = 1; itype <= ntypes; itype++) {
+      int n = strlen(arg[itype]) + 1;
+      typenames[itype] = new char[n];
+      strcpy(typenames[itype],arg[itype]);
+    }
+
+    return ntypes+1;
+  }
+
+  return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpXYZ::write_header(bigint n)
 {
   if (me == 0) {
     fprintf(fp,BIGINT_FORMAT "\n",n);
-    fprintf(fp,"Atoms\n");
+    fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-int DumpXYZ::count()
-{
-  if (igroup == 0) return atom->nlocal;
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  int m = 0;
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) m++;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void DumpXYZ::pack(int *ids)
 {
   int m,n;
@@ -112,7 +158,8 @@ void DumpXYZ::write_data(int n, double *mybuf)
   int m = 0;
   for (int i = 0; i < n; i++) {
     fprintf(fp,format,
-	    static_cast<int> (mybuf[m+1]),mybuf[m+2],mybuf[m+3],mybuf[m+4]);
+	    typenames[static_cast<int> (mybuf[m+1])],
+	    mybuf[m+2],mybuf[m+3],mybuf[m+4]);
     m += size_one;
   }
 }

src/dump_xyz.h

@@ -27,14 +27,17 @@ namespace LAMMPS_NS {
 class DumpXYZ : public Dump {
  public:
   DumpXYZ(class LAMMPS *, int, char**);
-  ~DumpXYZ() {}
+  ~DumpXYZ();
 
  private:
   void init_style();
   void write_header(bigint);
-  int count();
   void pack(int *);
   void write_data(int, double *);
+  int modify_param(int, char **);
+
+  int ntypes;
+  char **typenames;
 };
 
 }
@@ -55,4 +58,8 @@ E: Invalid dump xyz filename
 Filenames used with the dump xyz style cannot be binary or cause files
 to be written by each processor.
 
+E: Dump modify element names do not match atom types
+
+Number of element names must equal number of atom types.
+
 */