diff --git a/src/CLASS2/bond_class2.cpp b/src/CLASS2/bond_class2.cpp index 358375a6dc..097547c488 100644 --- a/src/CLASS2/bond_class2.cpp +++ b/src/CLASS2/bond_class2.cpp @@ -40,17 +40,15 @@ using int3_t = struct { int a, b, t; }; +using int4_t = struct { + int a, b, c, t; +}; + namespace LAMMPS_NS { -std::vector bond_class2_hot_i1; -std::vector bond_class2_hot_i2; -std::vector bond_class2_hot_delx; -std::vector bond_class2_hot_dely; -std::vector bond_class2_hot_delz; -std::vector bond_class2_hot_rsq; -std::vector bond_class2_hot_r; -bigint bond_class2_hot_timestep = -1; -int bond_class2_hot_nbondlist = 0; -bool bond_class2_hot_valid = false; +std::vector bond_class2_hot_angle_cache; +bigint bond_class2_hot_angle_timestep = -1; +int bond_class2_hot_angle_nanglelist = 0; +bool bond_class2_hot_angle_valid = false; } // namespace LAMMPS_NS /* ---------------------------------------------------------------------- */ @@ -83,7 +81,7 @@ void BondClass2::compute(int eflag, int vflag) int newton_bond = force->newton_bond; const int one_type = atom->nbondtypes == 1; - bond_class2_hot_valid = false; + bond_class2_hot_angle_valid = false; ev_init(eflag,vflag); @@ -128,20 +126,22 @@ void BondClass2::eval_one_type_hot(int nbondlist) const auto *_noalias const x = (dbl3_t *) atom->x[0]; auto *_noalias const f = (dbl3_t *) atom->f[0]; // NOLINT const int3_t *_noalias const bondlist = (int3_t *) neighbor->bondlist[0]; + const int nanglelist = neighbor->nanglelist; + const bool maybe_angle_cache = nbondlist == 2 * nanglelist; + const int4_t *_noalias const anglelist = maybe_angle_cache ? (int4_t *) neighbor->anglelist[0] : nullptr; const double coeff_r0 = r0[1]; const double coeff_2k2 = 2.0 * k2[1]; const double coeff_3k3 = 3.0 * k3[1]; const double coeff_4k4 = 4.0 * k4[1]; - bond_class2_hot_i1.resize(nbondlist); - bond_class2_hot_i2.resize(nbondlist); - bond_class2_hot_delx.resize(nbondlist); - bond_class2_hot_dely.resize(nbondlist); - bond_class2_hot_delz.resize(nbondlist); - bond_class2_hot_rsq.resize(nbondlist); - bond_class2_hot_r.resize(nbondlist); - bond_class2_hot_nbondlist = nbondlist; - bond_class2_hot_timestep = update->ntimestep; + bool angle_cache_valid = maybe_angle_cache; + if (maybe_angle_cache) { + bond_class2_hot_angle_cache.resize(nanglelist); + bond_class2_hot_angle_nanglelist = nanglelist; + bond_class2_hot_angle_timestep = update->ntimestep; + } else { + bond_class2_hot_angle_nanglelist = 0; + } int n = 0; @@ -164,20 +164,43 @@ void BondClass2::eval_one_type_hot(int nbondlist) const double r0 = sqrt(rsq0); const double r1 = sqrt(rsq1); - bond_class2_hot_i1[n] = i1_0; - bond_class2_hot_i2[n] = i2_0; - bond_class2_hot_delx[n] = delx0; - bond_class2_hot_dely[n] = dely0; - bond_class2_hot_delz[n] = delz0; - bond_class2_hot_rsq[n] = rsq0; - bond_class2_hot_r[n] = r0; - bond_class2_hot_i1[n + 1] = i1_1; - bond_class2_hot_i2[n + 1] = i2_1; - bond_class2_hot_delx[n + 1] = delx1; - bond_class2_hot_dely[n + 1] = dely1; - bond_class2_hot_delz[n + 1] = delz1; - bond_class2_hot_rsq[n + 1] = rsq1; - bond_class2_hot_r[n + 1] = r1; + if (angle_cache_valid) { + const int angle_index = n >> 1; + const int angle_i1 = anglelist[angle_index].a; + const int angle_i2 = anglelist[angle_index].b; + const int angle_i3 = anglelist[angle_index].c; + + if (i1_0 == angle_i2 && i1_1 == angle_i2) { + auto &cache = bond_class2_hot_angle_cache[angle_index]; + if (i2_0 == angle_i1 && i2_1 == angle_i3) { + cache.delx1 = -delx0; + cache.dely1 = -dely0; + cache.delz1 = -delz0; + cache.rsq1 = rsq0; + cache.r1 = r0; + cache.delx2 = -delx1; + cache.dely2 = -dely1; + cache.delz2 = -delz1; + cache.rsq2 = rsq1; + cache.r2 = r1; + } else if (i2_0 == angle_i3 && i2_1 == angle_i1) { + cache.delx1 = -delx1; + cache.dely1 = -dely1; + cache.delz1 = -delz1; + cache.rsq1 = rsq1; + cache.r1 = r1; + cache.delx2 = -delx0; + cache.dely2 = -dely0; + cache.delz2 = -delz0; + cache.rsq2 = rsq0; + cache.r2 = r0; + } else { + angle_cache_valid = false; + } + } else { + angle_cache_valid = false; + } + } const double dr0 = r0 - coeff_r0; const double dr1 = r1 - coeff_r0; @@ -226,14 +249,6 @@ void BondClass2::eval_one_type_hot(int nbondlist) const double rsq = delx * delx + dely * dely + delz * delz; const double r = sqrt(rsq); - bond_class2_hot_i1[n] = i1; - bond_class2_hot_i2[n] = i2; - bond_class2_hot_delx[n] = delx; - bond_class2_hot_dely[n] = dely; - bond_class2_hot_delz[n] = delz; - bond_class2_hot_rsq[n] = rsq; - bond_class2_hot_r[n] = r; - const double dr = r - coeff_r0; const double dr2 = dr * dr; const double dr3 = dr2 * dr; @@ -254,7 +269,7 @@ void BondClass2::eval_one_type_hot(int nbondlist) f[i2].z -= fz; } - bond_class2_hot_valid = true; + bond_class2_hot_angle_valid = angle_cache_valid; } template diff --git a/src/CLASS2/bond_class2.h b/src/CLASS2/bond_class2.h index d8ae006de7..6bcddc6e9e 100644 --- a/src/CLASS2/bond_class2.h +++ b/src/CLASS2/bond_class2.h @@ -22,8 +22,20 @@ BondStyle(class2,BondClass2); #include "bond.h" +#include + namespace LAMMPS_NS { +struct BondClass2AngleHotEntry { + double delx1, dely1, delz1, rsq1, r1; + double delx2, dely2, delz2, rsq2, r2; +}; + +extern std::vector bond_class2_hot_angle_cache; +extern bigint bond_class2_hot_angle_timestep; +extern int bond_class2_hot_angle_nanglelist; +extern bool bond_class2_hot_angle_valid; + class BondClass2 : public Bond { public: BondClass2(class LAMMPS *); diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp index af1e8d9451..0e42dc22eb 100644 --- a/src/MOLECULE/angle_harmonic.cpp +++ b/src/MOLECULE/angle_harmonic.cpp @@ -13,6 +13,7 @@ #include "angle_harmonic.h" +#include "../CLASS2/bond_class2.h" #include "atom.h" #include "comm.h" #include "domain.h" @@ -41,19 +42,6 @@ using int4_t = struct { int a, b, c, t; }; -namespace LAMMPS_NS { -extern std::vector bond_class2_hot_i1; -extern std::vector bond_class2_hot_i2; -extern std::vector bond_class2_hot_delx; -extern std::vector bond_class2_hot_dely; -extern std::vector bond_class2_hot_delz; -extern std::vector bond_class2_hot_rsq; -extern std::vector bond_class2_hot_r; -extern bigint bond_class2_hot_timestep; -extern int bond_class2_hot_nbondlist; -extern bool bond_class2_hot_valid; -} // namespace LAMMPS_NS - /* ---------------------------------------------------------------------- */ AngleHarmonic::AngleHarmonic(LAMMPS *_lmp) : Angle(_lmp) @@ -128,58 +116,119 @@ void AngleHarmonic::eval_one_type_hot(int nanglelist) const double coeff_k = k[1]; const double coeff_theta0 = theta0[1]; - if (bond_class2_hot_valid && bond_class2_hot_timestep == update->ntimestep && - bond_class2_hot_nbondlist == 2 * nanglelist) { - for (int n = 0; n < nanglelist; n++) { + if (bond_class2_hot_angle_valid && bond_class2_hot_angle_timestep == update->ntimestep && + bond_class2_hot_angle_nanglelist == nanglelist) { + int n = 0; + for (; n + 1 < nanglelist; n += 2) { + const int i1_0 = anglelist[n].a; + const int i2_0 = anglelist[n].b; + const int i3_0 = anglelist[n].c; + const int i1_1 = anglelist[n + 1].a; + const int i2_1 = anglelist[n + 1].b; + const int i3_1 = anglelist[n + 1].c; + const auto &_noalias cache0 = bond_class2_hot_angle_cache[n]; + const auto &_noalias cache1 = bond_class2_hot_angle_cache[n + 1]; + + const double delx1_0 = cache0.delx1; + const double dely1_0 = cache0.dely1; + const double delz1_0 = cache0.delz1; + const double delx2_0 = cache0.delx2; + const double dely2_0 = cache0.dely2; + const double delz2_0 = cache0.delz2; + const double rsq1_0 = cache0.rsq1; + const double rsq2_0 = cache0.rsq2; + const double r1_0 = cache0.r1; + const double r2_0 = cache0.r2; + + double c0 = delx1_0 * delx2_0 + dely1_0 * dely2_0 + delz1_0 * delz2_0; + c0 /= r1_0 * r2_0; + if (c0 > 1.0) c0 = 1.0; + if (c0 < -1.0) c0 = -1.0; + double s0 = sqrt(1.0 - c0 * c0); + if (s0 < SMALL) s0 = SMALL; + s0 = 1.0 / s0; + const double dtheta0 = acos(c0) - coeff_theta0; + const double tk0 = coeff_k * dtheta0; + const double a0 = -2.0 * tk0 * s0; + const double a11_0 = a0 * c0 / rsq1_0; + const double a12_0 = -a0 / (r1_0 * r2_0); + const double a22_0 = a0 * c0 / rsq2_0; + const double f1x0 = a11_0 * delx1_0 + a12_0 * delx2_0; + const double f1y0 = a11_0 * dely1_0 + a12_0 * dely2_0; + const double f1z0 = a11_0 * delz1_0 + a12_0 * delz2_0; + const double f3x0 = a22_0 * delx2_0 + a12_0 * delx1_0; + const double f3y0 = a22_0 * dely2_0 + a12_0 * dely1_0; + const double f3z0 = a22_0 * delz2_0 + a12_0 * delz1_0; + f[i1_0].x += f1x0; + f[i1_0].y += f1y0; + f[i1_0].z += f1z0; + f[i2_0].x -= f1x0 + f3x0; + f[i2_0].y -= f1y0 + f3y0; + f[i2_0].z -= f1z0 + f3z0; + f[i3_0].x += f3x0; + f[i3_0].y += f3y0; + f[i3_0].z += f3z0; + + const double delx1_1 = cache1.delx1; + const double dely1_1 = cache1.dely1; + const double delz1_1 = cache1.delz1; + const double delx2_1 = cache1.delx2; + const double dely2_1 = cache1.dely2; + const double delz2_1 = cache1.delz2; + const double rsq1_1 = cache1.rsq1; + const double rsq2_1 = cache1.rsq2; + const double r1_1 = cache1.r1; + const double r2_1 = cache1.r2; + + double c1 = delx1_1 * delx2_1 + dely1_1 * dely2_1 + delz1_1 * delz2_1; + c1 /= r1_1 * r2_1; + if (c1 > 1.0) c1 = 1.0; + if (c1 < -1.0) c1 = -1.0; + double s1 = sqrt(1.0 - c1 * c1); + if (s1 < SMALL) s1 = SMALL; + s1 = 1.0 / s1; + const double dtheta1 = acos(c1) - coeff_theta0; + const double tk1 = coeff_k * dtheta1; + const double a1 = -2.0 * tk1 * s1; + const double a11_1 = a1 * c1 / rsq1_1; + const double a12_1 = -a1 / (r1_1 * r2_1); + const double a22_1 = a1 * c1 / rsq2_1; + const double f1x1 = a11_1 * delx1_1 + a12_1 * delx2_1; + const double f1y1 = a11_1 * dely1_1 + a12_1 * dely2_1; + const double f1z1 = a11_1 * delz1_1 + a12_1 * delz2_1; + const double f3x1 = a22_1 * delx2_1 + a12_1 * delx1_1; + const double f3y1 = a22_1 * dely2_1 + a12_1 * dely1_1; + const double f3z1 = a22_1 * delz2_1 + a12_1 * delz1_1; + f[i1_1].x += f1x1; + f[i1_1].y += f1y1; + f[i1_1].z += f1z1; + f[i2_1].x -= f1x1 + f3x1; + f[i2_1].y -= f1y1 + f3y1; + f[i2_1].z -= f1z1 + f3z1; + f[i3_1].x += f3x1; + f[i3_1].y += f3y1; + f[i3_1].z += f3z1; + } + + for (; n < nanglelist; n++) { const int i1 = anglelist[n].a; const int i2 = anglelist[n].b; const int i3 = anglelist[n].c; - const int b0 = 2 * n; - const int b1 = b0 + 1; - - double delx1 = 0.0, dely1 = 0.0, delz1 = 0.0; - double delx2 = 0.0, dely2 = 0.0, delz2 = 0.0; - double rsq1 = 0.0, rsq2 = 0.0, r1 = 0.0, r2 = 0.0; - if (bond_class2_hot_i1[b0] == i2 && bond_class2_hot_i2[b0] == i1 && bond_class2_hot_i1[b1] == i2 && - bond_class2_hot_i2[b1] == i3) { - delx1 = -bond_class2_hot_delx[b0]; - dely1 = -bond_class2_hot_dely[b0]; - delz1 = -bond_class2_hot_delz[b0]; - rsq1 = bond_class2_hot_rsq[b0]; - r1 = bond_class2_hot_r[b0]; - delx2 = -bond_class2_hot_delx[b1]; - dely2 = -bond_class2_hot_dely[b1]; - delz2 = -bond_class2_hot_delz[b1]; - rsq2 = bond_class2_hot_rsq[b1]; - r2 = bond_class2_hot_r[b1]; - } else if (bond_class2_hot_i1[b0] == i2 && bond_class2_hot_i2[b0] == i3 && - bond_class2_hot_i1[b1] == i2 && bond_class2_hot_i2[b1] == i1) { - delx1 = -bond_class2_hot_delx[b1]; - dely1 = -bond_class2_hot_dely[b1]; - delz1 = -bond_class2_hot_delz[b1]; - rsq1 = bond_class2_hot_rsq[b1]; - r1 = bond_class2_hot_r[b1]; - delx2 = -bond_class2_hot_delx[b0]; - dely2 = -bond_class2_hot_dely[b0]; - delz2 = -bond_class2_hot_delz[b0]; - rsq2 = bond_class2_hot_rsq[b0]; - r2 = bond_class2_hot_r[b0]; - } else { - delx1 = x[i1].x - x[i2].x; - dely1 = x[i1].y - x[i2].y; - delz1 = x[i1].z - x[i2].z; - delx2 = x[i3].x - x[i2].x; - dely2 = x[i3].y - x[i2].y; - delz2 = x[i3].z - x[i2].z; - rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1; - rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2; - r1 = sqrt(rsq1); - r2 = sqrt(rsq2); - } + const auto &_noalias cache = bond_class2_hot_angle_cache[n]; + + const double delx1 = cache.delx1; + const double dely1 = cache.dely1; + const double delz1 = cache.delz1; + const double delx2 = cache.delx2; + const double dely2 = cache.dely2; + const double delz2 = cache.delz2; + const double rsq1 = cache.rsq1; + const double rsq2 = cache.rsq2; + const double r1 = cache.r1; + const double r2 = cache.r2; - const double r12 = r1 * r2; double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; - c /= r12; + c /= r1 * r2; if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; @@ -191,7 +240,7 @@ void AngleHarmonic::eval_one_type_hot(int nanglelist) const double tk = coeff_k * dtheta; const double a = -2.0 * tk * s; const double a11 = a * c / rsq1; - const double a12 = -a / r12; + const double a12 = -a / (r1 * r2); const double a22 = a * c / rsq2; const double f1x = a11 * delx1 + a12 * delx2;