diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp b/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp index 11d4aec318..2b86f87ec2 100644 --- a/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp +++ b/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) - ecoul = prefactor*erfc; + if (EFLAG) ecoul = prefactor*erfc; if (factor_coul < 1.0) { forcecoul -= (1.0-factor_coul)*prefactor; - if (EFLAG) - ecoul -= (1.0-factor_coul)*prefactor; + if (EFLAG) 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) - ecoul = qtmp*q[j] * table; + if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]); if (factor_coul < 1.0) { table = ctable[itable] + fraction*dctable[itable]; prefactor = qtmp*q[j] * table; forcecoul -= (1.0-factor_coul)*prefactor; - if (EFLAG) - ecoul -= (1.0-factor_coul)*prefactor; + if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor; } } } @@ -231,8 +227,7 @@ void PairLJSDKCoulLong::eval() - lj4[itype][jtype]) - offset[itype][jtype]; } forcelj *= factor_lj; - if (EFLAG) - evdwl *= factor_lj; + if (EFLAG) evdwl *= factor_lj; } fpair = (forcecoul + forcelj) * r2inv; @@ -247,8 +242,7 @@ void PairLJSDKCoulLong::eval() } if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR, - evdwl,ecoul,fpair,delx,dely,delz); - + evdwl,ecoul,fpair,delx,dely,delz); } } f[i][0] += fxtmp; diff --git a/src/USER-OMP/improper_class2_omp.cpp b/src/USER-OMP/improper_class2_omp.cpp index b59cf9e4b1..4204502396 100644 --- a/src/USER-OMP/improper_class2_omp.cpp +++ b/src/USER-OMP/improper_class2_omp.cpp @@ -79,7 +79,7 @@ void ImproperClass2OMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void ImproperClass2OMP::eval(int nfrom, int nto, ThrData * const thr) { int i1,i2,i3,i4,i,j,k,n,type; diff --git a/src/USER-OMP/pair_buck_coul_omp.cpp b/src/USER-OMP/pair_buck_coul_omp.cpp index 7a66fa2fb0..4ff10ba02c 100644 --- a/src/USER-OMP/pair_buck_coul_omp.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_coul_long_omp.cpp b/src/USER-OMP/pair_coul_long_omp.cpp index 62a76783a7..8be5e35a0f 100644 --- a/src/USER-OMP/pair_coul_long_omp.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_lj96_cut_omp.cpp b/src/USER-OMP/pair_lj96_cut_omp.cpp index 69d9539acf..b92749bf56 100644 --- a/src/USER-OMP/pair_lj96_cut_omp.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_lj_charmm_coul_long_omp.cpp b/src/USER-OMP/pair_lj_charmm_coul_long_omp.cpp index 680626daa5..9a5b1300e3 100644 --- a/src/USER-OMP/pair_lj_charmm_coul_long_omp.cpp +++ b/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 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; - numneigh = list->numneigh; - firstneigh = list->firstneigh; + 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]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + const int i = ilist[ii]; + const int itype = type[i]; + const double qtmp = q[i]; + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + double fxtmp,fytmp,fztmp; fxtmp=fytmp=fztmp=0.0; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; + const int * const jlist = firstneigh[i]; + const int jnum = numneigh[i]; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; + for (int jj = 0; jj < jnum; jj++) { + double forcecoul, forcelj, evdwl, ecoul; + forcecoul = forcelj = evdwl = ecoul = 0.0; + + 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); - 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; + 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 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; - itable >>= ncoulshiftbits; - fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; + const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits; + const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + const double 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; + if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]); + if (sbindex) { + const double table2 = ctable[itable] + fraction*dctable[itable]; + 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; - jtype = type[j]; + const double r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[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; - switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; - philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); - forcelj = forcelj*switch1 + philj*switch2; - } - forcelj *= factor_lj; - } else forcelj = 0.0; + if (EFLAG) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); - fpair = (forcecoul + forcelj) * r2inv; + if (rsq > cut_lj_innersq) { + const double drsq = cut_ljsq - rsq; + const double cut2 = (rsq - cut_lj_innersq) * drsq; + const double switch1 = drsq * (drsq*drsq + 3.0*cut2) * inv_denom_lj; + const double switch2 = 12.0*rsq * cut2 * inv_denom_lj; + 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; + } + } + + 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 (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, + if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR, evdwl,ecoul,fpair,delx,dely,delz,thr); } } diff --git a/src/USER-OMP/pair_lj_coul_omp.cpp b/src/USER-OMP/pair_lj_coul_omp.cpp index 6f12ab2d84..28e4a47645 100644 --- a/src/USER-OMP/pair_lj_coul_omp.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_lj_cut_coul_long_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_long_omp.cpp index 238dc6f386..5f9b80ac80 100644 --- a/src/USER-OMP/pair_lj_cut_coul_long_omp.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp index b9d83250c3..3e03fa6b97 100644 --- a/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp @@ -145,7 +145,7 @@ void PairLJCutCoulLongTIP4POMP::compute(int eflag, int vflag) /* ---------------------------------------------------------------------- */ -template +template void PairLJCutCoulLongTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,ii,jj,jnum,itype,jtype,itable; diff --git a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp index d4c0362723..eb780621cb 100644 --- a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp @@ -167,7 +167,7 @@ void PairLJCutCoulPPPMTIP4POMP::compute(int eflag, int vflag) /* ---------------------------------------------------------------------- */ -template +template void PairLJCutCoulPPPMTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,ii,jj,jnum,itype,jtype,itable; diff --git a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.h b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.h index f4a0bf1f53..4596369324 100644 --- a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.h +++ b/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 + template void eval(int ifrom, int ito, ThrData * const thr); class PPPMTIP4PProxy *kspace; diff --git a/src/USER-OMP/pair_lj_cut_omp.cpp b/src/USER-OMP/pair_lj_cut_omp.cpp index 0c86e1be9b..f8d3b1f8ac 100644 --- a/src/USER-OMP/pair_lj_cut_omp.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_lj_sdk_coul_long_omp.cpp b/src/USER-OMP/pair_lj_sdk_coul_long_omp.cpp index 580c57d96c..6632f4fc47 100644 --- a/src/USER-OMP/pair_lj_sdk_coul_long_omp.cpp +++ b/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 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]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; + const int i = ilist[ii]; + const int itype = type[i]; + const double qtmp = q[i]; + const double xtmp = x[i][0]; + 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]; - factor_lj = special_lj[sbmask(j)]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; + const int sbindex = sbmask(jlist[jj]); + const int j = jlist[jj] & NEIGHMASK; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; + 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; const int ljt = lj_type[itype][jtype]; 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; - prefactor = qqrd2e * qtmp*q[j]/r; + 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 prefactor = qqrd2e * qtmp*q[j]/r; forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (EFLAG) - ecoul = prefactor*erfc; - if (factor_coul < 1.0) { - forcecoul -= (1.0-factor_coul)*prefactor; - if (EFLAG) - ecoul -= (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; - itable >>= ncoulshiftbits; - fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; + const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits; + const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + const double table = ftable[itable] + fraction*dftable[itable]; forcecoul = qtmp*q[j] * table; - if (EFLAG) - 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) - ecoul -= (1.0-factor_coul)*prefactor; + if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]); + if (sbindex) { + const double table2 = ctable[itable] + fraction*dctable[itable]; + const double prefactor = qtmp*q[j] * table2; + const double adjust = (1.0-special_coul[sbindex])*prefactor; + forcecoul -= adjust; + if (EFLAG) ecoul -= adjust; } } } @@ -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]; } - forcelj *= factor_lj; - if (EFLAG) - evdwl *= factor_lj; + + if (sbindex) { + const double factor_lj = special_lj[sbindex]; + forcelj *= factor_lj; + if (EFLAG) 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; diff --git a/src/XTC/dump_xtc.cpp b/src/XTC/dump_xtc.cpp index e60f5e0670..d003a1c20f 100644 --- a/src/XTC/dump_xtc.cpp +++ b/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; diff --git a/src/XTC/dump_xtc.h b/src/XTC/dump_xtc.h index d6aee90bec..e0a113e09d 100644 --- a/src/XTC/dump_xtc.h +++ b/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(); diff --git a/src/dump.cpp b/src/dump.cpp index feda8ea156..f5dcea53c3 100644 --- a/src/dump.cpp +++ b/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 diff --git a/src/dump.h b/src/dump.h index e6b5460287..852202544a 100644 --- a/src/dump.h +++ b/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; diff --git a/src/dump_atom.cpp b/src/dump_atom.cpp index a9dfdcc03c..13930a50c6 100644 --- a/src/dump_atom.cpp +++ b/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); diff --git a/src/dump_atom.h b/src/dump_atom.h index 682b660049..54a57784dc 100644 --- a/src/dump_atom.h +++ b/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 *); diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp index 484c8056f1..6512229735 100644 --- a/src/dump_dcd.cpp +++ b/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; diff --git a/src/dump_dcd.h b/src/dump_dcd.h index a0197021b9..6915068579 100644 --- a/src/dump_dcd.h +++ b/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 **); diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index 817e376fbf..8d9cd6b151 100644 --- a/src/dump_xyz.cpp +++ b/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 (mybuf[m+1]),mybuf[m+2],mybuf[m+3],mybuf[m+4]); + typenames[static_cast (mybuf[m+1])], + mybuf[m+2],mybuf[m+3],mybuf[m+4]); m += size_one; } } diff --git a/src/dump_xyz.h b/src/dump_xyz.h index 57016eafc8..4824ed6df4 100644 --- a/src/dump_xyz.h +++ b/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. + */