From d02ff1d70efce95f9487d69f025e2e56092d431d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 19 Aug 2022 16:33:55 -0400 Subject: [PATCH] correct force and energy tally for using full neighbor lists and newton on/off --- src/DIELECTRIC/pair_coul_cut_dielectric.cpp | 8 ---- src/DIELECTRIC/pair_coul_long_dielectric.cpp | 8 ---- .../pair_lj_cut_coul_cut_dielectric.cpp | 17 +++------ .../pair_lj_cut_coul_debye_dielectric.cpp | 8 ---- .../pair_lj_cut_coul_long_dielectric.cpp | 9 ----- .../pair_lj_cut_coul_msm_dielectric.cpp | 18 --------- .../pair_lj_long_coul_long_dielectric.cpp | 9 ----- .../pair_lj_cut_coul_long_dielectric_omp.cpp | 6 --- src/OPENMP/thr_omp.cpp | 37 +++++++++++++++++++ src/OPENMP/thr_omp.h | 2 + 10 files changed, 44 insertions(+), 78 deletions(-) diff --git a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp index aaa39f6aea..7c6c4cd11e 100644 --- a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp @@ -137,16 +137,8 @@ void PairCoulCutDielectric::compute(int eflag, int vflag) efield[i][1] += dely * efield_i; efield[i][2] += delz * efield_i; - if (newton_pair && j >= nlocal) { - fpair_j = factor_coul * eps[j] * forcecoul * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } - if (eflag) { ecoul = factor_coul * qqrd2e * scale[itype][jtype] * qtmp * q[j] * (etmp + eps[j]) * rinv; - ecoul *= 0.5; } if (evflag) ev_tally_full(i, 0.0, ecoul, fpair_i, delx, dely, delz); } diff --git a/src/DIELECTRIC/pair_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_coul_long_dielectric.cpp index 620f0cca34..ddfff752f4 100644 --- a/src/DIELECTRIC/pair_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_long_dielectric.cpp @@ -169,13 +169,6 @@ void PairCoulLongDielectric::compute(int eflag, int vflag) efield[i][1] += dely * efield_i; efield[i][2] += delz * efield_i; - if (newton_pair && j >= nlocal) { - fpair_j = eps[j] * forcecoul * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } - if (eflag) { if (!ncoultablebits || rsq <= tabinnersq) ecoul = prefactor * (etmp + eps[j]) * erfc; @@ -183,7 +176,6 @@ void PairCoulLongDielectric::compute(int eflag, int vflag) table = etable[itable] + fraction * detable[itable]; ecoul = scale[itype][jtype] * qtmp * q[j] * (etmp + eps[j]) * table; } - ecoul *= 0.5; if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; } diff --git a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp index 517fbc481a..243bf5b2e9 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp @@ -57,7 +57,7 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag) int i, j, ii, jj, inum, jnum, itype, jtype; double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul; double fpair_i, fpair_j; - double rsq, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj, efield_i, epot_i; + double rsq, rinv, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj, efield_i, epot_i; int *ilist, *jlist, *numneigh, **firstneigh; if (atom->nmax > nmax) { @@ -114,7 +114,7 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag) efield[i][0] = efield[i][1] = efield[i][2] = 0; } - epot[i] = 0; + epot[i] = 0.0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -130,9 +130,10 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag) if (rsq < cutsq[itype][jtype]) { r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); if (rsq < cut_coulsq[itype][jtype] && rsq > EPSILON) { - efield_i = qqrd2e * q[j] * sqrt(r2inv); + efield_i = qqrd2e * q[j] * rinv; forcecoul = qtmp * efield_i; epot_i = efield_i; } else @@ -156,19 +157,11 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag) epot[i] += epot_i; - if (newton_pair && j >= nlocal) { - fpair_j = (factor_coul * eps[j] * forcecoul + factor_lj * forcelj) * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } - if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { - ecoul = factor_coul * qqrd2e * qtmp * q[j] * (etmp + eps[j]) * sqrt(r2inv); + ecoul = factor_coul * qqrd2e * qtmp * q[j] * (etmp + eps[j]) *rinv; } else ecoul = 0.0; - ecoul *= 0.5; if (rsq < cut_ljsq[itype][jtype]) { evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; diff --git a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp index 814215872b..53de610b9a 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp @@ -159,19 +159,11 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag) efield[i][2] += delz * efield_i; epot[i] += epot_i; - if (newton_pair && j >= nlocal) { - fpair_j = (factor_coul * eps[j] * forcecoul + factor_lj * forcelj) * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } - if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { ecoul = factor_coul * qqrd2e * qtmp * q[j] * (etmp + eps[j]) * rinv * screening; } else ecoul = 0.0; - ecoul *= 0.5; if (rsq < cut_ljsq[itype][jtype]) { evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; diff --git a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp index dad8d74617..8c6b854b1d 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp @@ -196,14 +196,6 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag) epot[i] += epot_i; - if (newton_pair && j >= nlocal) { - - fpair_j = (forcecoul * eps[j] + factor_lj * forcelj) * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } - if (eflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) @@ -212,7 +204,6 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag) table = etable[itable] + fraction * detable[itable]; ecoul = qtmp * q[j] * (etmp + eps[j]) * table; } - ecoul *= 0.5; if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; } else ecoul = 0.0; diff --git a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp index 0c66f52a4f..6fbd8c400c 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp @@ -208,12 +208,6 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag) efield[i][1] += dely * efield_i; efield[i][2] += delz * efield_i; - if (newton_pair && j >= nlocal) { - fpair_j = (forcecoul * eps[j] + factor_lj * forcelj) * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } } else { // separate LJ and Coulombic forces @@ -223,11 +217,6 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag) f[i][0] += delx * fpair; f[i][1] += dely * fpair; f[i][2] += delz * fpair; - if (newton_pair) { - f[j][0] -= delx * fpair; - f[j][1] -= dely * fpair; - f[j][2] -= delz * fpair; - } fpair_i = (forcecoul * etmp) * r2inv; ftmp[i][0] += delx * fpair_i; @@ -238,13 +227,6 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag) efield[i][0] += delx * efield_i; efield[i][1] += dely * efield_i; efield[i][2] += delz * efield_i; - - if (newton_pair && j >= nlocal) { - fpair_j = (forcecoul * eps[j]) * r2inv; - ftmp[j][0] -= delx * fpair_j; - ftmp[j][1] -= dely * fpair_j; - ftmp[j][2] -= delz * fpair_j; - } } if (eflag) { diff --git a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp index f99f6438ea..0ceacb206d 100644 --- a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp @@ -268,14 +268,6 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag) epot[i] += epot_i; - if (newton_pair && j >= nlocal) { - - fpair_j = (force_coul * eps[j] + factor_lj * force_lj) * r2inv; - f[j][0] -= delx * fpair_j; - f[j][1] -= dely * fpair_j; - f[j][2] -= delz * fpair_j; - } - if (eflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) @@ -284,7 +276,6 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag) table = etable[itable] + fraction * detable[itable]; ecoul = qtmp * q[j] * (etmp + eps[j]) * table; } - ecoul *= 0.5; if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; } else ecoul = 0.0; diff --git a/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp index 8e3880810c..f197a9a3b4 100644 --- a/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp +++ b/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp @@ -224,12 +224,6 @@ void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData *const t eztmp += delz * efield_i; epot[i] += epot_i; - if (NEWTON_PAIR || j < nlocal) { - f[j].x -= delx * fpair; - f[j].y -= dely * fpair; - f[j].z -= delz * fpair; - } - if (EFLAG) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) diff --git a/src/OPENMP/thr_omp.cpp b/src/OPENMP/thr_omp.cpp index 95392b17a4..cc7efdf603 100644 --- a/src/OPENMP/thr_omp.cpp +++ b/src/OPENMP/thr_omp.cpp @@ -588,6 +588,43 @@ void ThrOMP::ev_tally_thr(Pair * const pair, const int i, const int j, const int } } +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into per thread global and per-atom accumulators + for full neighbor list. Only tally atom i and also set newton to off. +------------------------------------------------------------------------- */ + +void ThrOMP::ev_tally_full_thr(Pair * const pair, const int i, const double evdwl, + const double ecoul, const double fpair, const double delx, + const double dely, const double delz, ThrData * const thr) +{ + + if (pair->eflag_either) + e_tally_thr(pair, i, i+1, i, 0, evdwl, ecoul, thr); + + if (pair->vflag_either) { + double v[6]; + v[0] = delx*delx*fpair; + v[1] = dely*dely*fpair; + v[2] = delz*delz*fpair; + v[3] = delx*dely*fpair; + v[4] = delx*delz*fpair; + v[5] = dely*delz*fpair; + + v_tally_thr(pair, i, i+1, i, 0, v, thr); + } + + if (pair->num_tally_compute > 0) { + // ev_tally callbacks are not thread safe and thus have to be protected +#if defined(_OPENMP) +#pragma omp critical +#endif + for (int k=0; k < pair->num_tally_compute; ++k) { + Compute *c = pair->list_tally_compute[k]; + c->pair_tally_callback(i, i+1, i, 0, evdwl, ecoul, fpair, delx, dely, delz); + } + } +} + /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz diff --git a/src/OPENMP/thr_omp.h b/src/OPENMP/thr_omp.h index 78d93bd55f..b1ecf93a88 100644 --- a/src/OPENMP/thr_omp.h +++ b/src/OPENMP/thr_omp.h @@ -129,6 +129,8 @@ class ThrOMP { void ev_tally_thr(Pair *const, const int, const int, const int, const int, const double, const double, const double, const double, const double, const double, ThrData *const); + void ev_tally_full_thr(Pair *const, const int, const double, const double, const double, + const double, const double, const double, ThrData *const); void ev_tally_xyz_thr(Pair *const, const int, const int, const int, const int, const double, const double, const double, const double, const double, const double, const double, const double, ThrData *const);