From 85765a2bf3242a2f6deb14a2d8c7084ec13076b4 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 12 Jun 2023 12:45:52 +0300 Subject: [PATCH 01/60] Include born_matrix() definition in bond_gaussian.h --- src/EXTRA-MOLECULE/bond_gaussian.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/bond_gaussian.h b/src/EXTRA-MOLECULE/bond_gaussian.h index 7af6f1f4d9..e466df47d4 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.h +++ b/src/EXTRA-MOLECULE/bond_gaussian.h @@ -35,6 +35,7 @@ class BondGaussian : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; protected: int *nterms; From a05fcc326efd1f7ebef08516ddf15b262c155215 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 12 Jun 2023 12:47:21 +0300 Subject: [PATCH 02/60] Implement born_matrix() in bond_gaussian.cpp --- src/EXTRA-MOLECULE/bond_gaussian.cpp | 45 +++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-MOLECULE/bond_gaussian.cpp b/src/EXTRA-MOLECULE/bond_gaussian.cpp index baca0b6e1a..816ab51516 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.cpp +++ b/src/EXTRA-MOLECULE/bond_gaussian.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ #include "bond_gaussian.h" - +#include #include "atom.h" #include "comm.h" #include "error.h" @@ -35,6 +35,7 @@ BondGaussian::BondGaussian(LAMMPS *lmp) : Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr), r0(nullptr) { + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -294,3 +295,45 @@ double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/, double & return -(force->boltz * bond_temperature[type]) * log(sum_g_i); } + +/* ---------------------------------------------------------------------- */ + +void BondGaussian::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + + // first derivative of energy with respect to distance + double sum_g_i = 0.0; + double sum_numerator = 0.0; + for (int i = 0; i < nterms[type]; i++) { + double dr = r - r0[type][i]; + double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2))); + double exponent = -2 * dr * dr / (width[type][i] * width[type][i]); + double g_i = prefactor * exp(exponent); + sum_g_i += g_i; + sum_numerator += g_i * dr / (width[type][i] * width[type][i]); + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + du = 4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i); + + // second derivative of energy with respect to distance + sum_g_i = 0.0; + double sum_dg_i = 0.0; + double sum_d2g_i = 0.0; + for (int i = 0; i < nterms[type]; i++) { + double dr = r - r0[type][i]; + double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2))); + double exponent = -2 * dr * dr / (width[type][i] * width[type][i]); + double g_i = prefactor * exp(exponent); + sum_g_i += g_i; + sum_dg_i -= 4.0 * g_i * dr / pow(width[type][i], 2); + sum_d2g_i += 4.0 * g_i * (4.0 * pow(r0[type][i], 2) - 8.0 * r0[type][i] * r - pow(width[type][i], 2) + 4.0 * r * r) / pow(width[type][i], 4) ; + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + double numerator = sum_d2g_i*sum_g_i - sum_dg_i*sum_dg_i; + double denominator = sum_g_i * sum_g_i; + + du2 = - (force->boltz * bond_temperature[type]) * numerator / denominator; +} From 42c843ff4f783069d409404be42fa89c43f60ccc Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 18 Jun 2023 18:24:24 +0300 Subject: [PATCH 03/60] remove iostream from bond_gaussian.cpp --- src/EXTRA-MOLECULE/bond_gaussian.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-MOLECULE/bond_gaussian.cpp b/src/EXTRA-MOLECULE/bond_gaussian.cpp index 816ab51516..9a8546e278 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.cpp +++ b/src/EXTRA-MOLECULE/bond_gaussian.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ #include "bond_gaussian.h" -#include + #include "atom.h" #include "comm.h" #include "error.h" From 4e17cc551e93faf1433bb002bfbde9838684c2eb Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 18 Jun 2023 18:25:35 +0300 Subject: [PATCH 04/60] inlcude born_matrix() definition in in angle_quartic.h --- src/EXTRA-MOLECULE/angle_quartic.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/angle_quartic.h b/src/EXTRA-MOLECULE/angle_quartic.h index 3f0396f27b..7de51b24d1 100644 --- a/src/EXTRA-MOLECULE/angle_quartic.h +++ b/src/EXTRA-MOLECULE/angle_quartic.h @@ -35,6 +35,7 @@ class AngleQuartic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; protected: double *k2, *k3, *k4, *theta0; From eb8512ba2a988b9baea384fe96aa935c28e6005a Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 18 Jun 2023 18:26:48 +0300 Subject: [PATCH 05/60] implementation of born_matrix() for angle_quartic.cpp --- src/EXTRA-MOLECULE/angle_quartic.cpp | 41 +++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-MOLECULE/angle_quartic.cpp b/src/EXTRA-MOLECULE/angle_quartic.cpp index f28e209a77..eaccdbe608 100644 --- a/src/EXTRA-MOLECULE/angle_quartic.cpp +++ b/src/EXTRA-MOLECULE/angle_quartic.cpp @@ -37,7 +37,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -AngleQuartic::AngleQuartic(LAMMPS *lmp) : Angle(lmp) {} +AngleQuartic::AngleQuartic(LAMMPS *lmp) : Angle(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -286,3 +289,39 @@ double AngleQuartic::single(int type, int i1, int i2, int i3) double dtheta4 = dtheta3 * dtheta; return k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4; } + +/* ---------------------------------------------------------------------- */ + +void AngleQuartic::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c); + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + double dtheta = theta - theta0[type]; + double dtheta2 = dtheta * dtheta; + double dtheta3 = dtheta2 * dtheta; + + du = -(2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3) / s; + du2 = (2.0 * k2[type] + 6.0 * k3[type] * dtheta + 12.0 * k4[type] * dtheta2) / (s*s) - + (2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3) * c / (s*s*s); +} From 365f4bc55945506de57e44dc0c6ec7102fa2d000 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 18 Jun 2023 18:44:28 +0300 Subject: [PATCH 06/60] non-zero born_matrix_enable flag in angle_fourier.cpp --- src/EXTRA-MOLECULE/angle_fourier.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/angle_fourier.cpp b/src/EXTRA-MOLECULE/angle_fourier.cpp index 549da0c196..c7eb3d4fe4 100644 --- a/src/EXTRA-MOLECULE/angle_fourier.cpp +++ b/src/EXTRA-MOLECULE/angle_fourier.cpp @@ -39,6 +39,7 @@ using namespace MathConst; AngleFourier::AngleFourier(LAMMPS *lmp) : Angle(lmp) { + born_matrix_enable = 1; k = nullptr; C0 = nullptr; C1 = nullptr; From 2f227614615219d13538481646563c2a8111f251 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 19 Jun 2023 16:38:34 +0300 Subject: [PATCH 07/60] born_matrix() method in bond_mm3.h --- src/YAFF/bond_mm3.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/YAFF/bond_mm3.h b/src/YAFF/bond_mm3.h index 302c4052d0..ea89ac826d 100644 --- a/src/YAFF/bond_mm3.h +++ b/src/YAFF/bond_mm3.h @@ -35,6 +35,7 @@ class BondMM3 : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; protected: double *r0, *k2; From bfc969d5c5ee4da4e9062bffe7f96d55ffb318ab Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 19 Jun 2023 16:39:49 +0300 Subject: [PATCH 08/60] implementation of born_matrix in bond_mm3.cpp --- src/YAFF/bond_mm3.cpp | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/YAFF/bond_mm3.cpp b/src/YAFF/bond_mm3.cpp index a5ef6fb8bc..31ce2dad3e 100644 --- a/src/YAFF/bond_mm3.cpp +++ b/src/YAFF/bond_mm3.cpp @@ -31,7 +31,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondMM3::BondMM3(LAMMPS *lmp) : Bond(lmp) {} +BondMM3::BondMM3(LAMMPS *lmp) : Bond(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -219,3 +222,19 @@ double BondMM3::single(int type, double rsq, else fforce = 0.0; return k2[type]*dr2*(1.0+K3*dr+K4*dr2); } + +/* ---------------------------------------------------------------------- */ + +void BondMM3::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + double dr = r - r0[type]; + double dr2 = dr * dr; + double dr3 = dr2 * dr; + + double K3 = -2.55 * k2[type] /force->angstrom; + double K4 = 7.0 * k2[type] * 2.55 * 2.55 / (12.0 * force->angstrom * force->angstrom); + + du = 2.0 * k2[type] * dr + 3.0 * K3 * dr2 + 4.0 * K4 * dr3; + du2 = 2.0 * k2[type] + 6.0 * K3 * dr + 12.0 * K4 * dr2; +} From 6c9d42b7c363a66fa56739f7d3ba3dbc5a1952c2 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 19 Jun 2023 19:38:50 +0300 Subject: [PATCH 09/60] Include born_matrix() definition in bond_harmonic_shift_cut.h --- src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h index 752ac010d9..09d6ab5330 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h +++ b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h @@ -35,6 +35,7 @@ class BondHarmonicShiftCut : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; protected: double *k, *r0, *r1; From f6b259b1866d9d5316d6464eae7d551fea9fc760 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 19 Jun 2023 19:40:11 +0300 Subject: [PATCH 10/60] Implementing born_matrix in bond_harmonic_shift_cut.cpp --- .../bond_harmonic_shift_cut.cpp | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp index fedcb95ee8..ebcfdb0258 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp @@ -31,7 +31,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondHarmonicShiftCut::BondHarmonicShiftCut(LAMMPS *lmp) : Bond(lmp) {} +BondHarmonicShiftCut::BondHarmonicShiftCut(LAMMPS *lmp) : Bond(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -219,3 +222,19 @@ double BondHarmonicShiftCut::single(int type, double rsq, int /*i*/, int /*j*/, fforce = -2.0*k[type]*dr/r; return k[type]*(dr*dr - dr2*dr2); } + +/* ---------------------------------------------------------------------- */ + +void BondHarmonicShiftCut::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + du = 0.0; + du2 = 0.0; + + double r = sqrt(rsq); + if (r>r1[type]) return; + + double dr = r - r0[type]; + + du2 = 2 * k[type]; + if (r > 0.0) du = du2 * dr; +} From ad3752431fe9711ef40f2d53cd28d7ab83153b73 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 19 Jun 2023 19:42:01 +0300 Subject: [PATCH 11/60] Regular pointer for coord and coordp in compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index 2b0ffef0f8..1b5b3a911f 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -49,7 +49,7 @@ class ComputeStressMopProfile : public Compute { int originflag; double origin, delta, offset, invdelta; int nbins; - double **coord, **coordp; + double *coord, *coordp; double **values_local, **values_global; double **bond_local, **bond_global; double **local_contribution; From dc1eb43cf2f15c61c0b44bb872407c4b5091e135 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 19 Jun 2023 19:47:34 +0300 Subject: [PATCH 12/60] Cleaning coord and coordp vectors in compute_stress_mop_profile.cpp --- .../compute_stress_mop_profile.cpp | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 1f1dc30733..3a85869f3c 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -258,7 +258,7 @@ void ComputeStressMopProfile::compute_array() MPI_Allreduce(&bond_local[0][0],&bond_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); for (int ibin=0; ibin ((hi-lo) * invdelta + 1.5); //allocate bin arrays - memory->create(coord,nbins,1,"stress/mop/profile:coord"); - memory->create(coordp,nbins,1,"stress/mop/profile:coordp"); + memory->create(coord,nbins,"stress/mop/profile:coord"); + memory->create(coordp,nbins,"stress/mop/profile:coordp"); memory->create(values_local,nbins,nvalues,"stress/mop/profile:values_local"); memory->create(values_global,nbins,nvalues,"stress/mop/profile:values_global"); memory->create(bond_local,nbins,nvalues,"stress/mop/profile:bond_local"); @@ -652,11 +652,11 @@ void ComputeStressMopProfile::setup_bins() // set bin coordinates for (i = 0; i < nbins; i++) { - coord[i][0] = offset + i*delta; - if (coord[i][0] < (domain->boxlo[dir]+domain->prd_half[dir])) { - coordp[i][0] = coord[i][0] + domain->prd[dir]; + coord[i] = offset + i*delta; + if (coord[i] < (domain->boxlo[dir]+domain->prd_half[dir])) { + coordp[i] = coord[i] + domain->prd[dir]; } else { - coordp[i][0] = coord[i][0] - domain->prd[dir]; + coordp[i] = coord[i] - domain->prd[dir]; } } } From 2631a159affeedef3dd2eeaaa1c3ca0529396b97 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 20 Jun 2023 15:41:09 +0300 Subject: [PATCH 13/60] define born_matrix in angle_mm3.h --- src/YAFF/angle_mm3.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/YAFF/angle_mm3.h b/src/YAFF/angle_mm3.h index 95009a9cf6..22f5bd746c 100644 --- a/src/YAFF/angle_mm3.h +++ b/src/YAFF/angle_mm3.h @@ -35,6 +35,7 @@ class AngleMM3 : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; protected: double *theta0, *k2; From bb2d691e7863e5b4746da6fd8f6a951581cbde51 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 20 Jun 2023 15:42:47 +0300 Subject: [PATCH 14/60] implement born_matrix in angle_mm3.cpp --- src/YAFF/angle_mm3.cpp | 45 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/src/YAFF/angle_mm3.cpp b/src/YAFF/angle_mm3.cpp index c75a0d8308..af199f6fe9 100644 --- a/src/YAFF/angle_mm3.cpp +++ b/src/YAFF/angle_mm3.cpp @@ -36,7 +36,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -AngleMM3::AngleMM3(LAMMPS *lmp) : Angle(lmp) {} +AngleMM3::AngleMM3(LAMMPS *lmp) : Angle(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -284,3 +287,43 @@ double AngleMM3::single(int type, int i1, int i2, int i3) return energy; } + +/* ---------------------------------------------------------------------- */ + +void AngleMM3::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c); + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + double dtheta = theta - theta0[type]; + double dtheta2 = dtheta*dtheta; + double dtheta3 = dtheta2*dtheta; + double dtheta4 = dtheta3*dtheta; + double dtheta5 = dtheta4*dtheta; + double df = 2.0 * dtheta - 2.406423 * dtheta2 + 0.735348 * dtheta3 - 0.65832 * dtheta4 + 1.42254 * dtheta5; + double d2f = 2.0 - 4.812846 * dtheta + 2.206044 * dtheta2 - 2.63328 * dtheta3 + 7.1127 * dtheta4; + + du = -k2[type] * df / s; + du2 = k2[type] * (d2f - df * c / s) / (s * s) ; +} From 345a834c7e0693770e37a13005118ad29e19ff7b Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 20 Jun 2023 16:20:57 +0300 Subject: [PATCH 15/60] Include definition of born_matrix() in angle_cosine_periodic.h --- src/EXTRA-MOLECULE/angle_cosine_periodic.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/angle_cosine_periodic.h b/src/EXTRA-MOLECULE/angle_cosine_periodic.h index 4e584b4543..f04ed04784 100644 --- a/src/EXTRA-MOLECULE/angle_cosine_periodic.h +++ b/src/EXTRA-MOLECULE/angle_cosine_periodic.h @@ -35,6 +35,7 @@ class AngleCosinePeriodic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; protected: double *k; From 7f3a930d8923a13fb14c02e85d56edb703b1c2cd Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 20 Jun 2023 16:21:57 +0300 Subject: [PATCH 16/60] Implement born_matrix() in angle_cosine_periodic.cpp --- src/EXTRA-MOLECULE/angle_cosine_periodic.cpp | 40 +++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp b/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp index 15d0575f6d..34a8e9d8e5 100644 --- a/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp +++ b/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp @@ -38,7 +38,10 @@ using namespace MathSpecial; /* ---------------------------------------------------------------------- */ -AngleCosinePeriodic::AngleCosinePeriodic(LAMMPS *lmp) : Angle(lmp) {} +AngleCosinePeriodic::AngleCosinePeriodic(LAMMPS *lmp) : Angle(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -298,3 +301,38 @@ double AngleCosinePeriodic::single(int type, int i1, int i2, int i3) c = cos(acos(c)*multiplicity[type]); return 2.0*k[type]*(1.0-b[type]*powsign(multiplicity[type])*c); } + +/* ---------------------------------------------------------------------- */ + +void AngleCosinePeriodic::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c); + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + double m_angle = multiplicity[type] * theta; + double prefactor = -2.0 * k[type] * b[type] * powsign(multiplicity[type]) * multiplicity[type]; + + du = prefactor * sin(m_angle) / s; + du2 = prefactor * (c * sin(m_angle) - s * cos(m_angle) * multiplicity[type]) / (s * s * s); +} From fb31ffe17cbff6854bc60870a8be54258df7ed07 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 21 Jun 2023 10:56:53 +0300 Subject: [PATCH 17/60] Definition of born_matrix() in dihedral_helix.h --- src/EXTRA-MOLECULE/dihedral_helix.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/dihedral_helix.h b/src/EXTRA-MOLECULE/dihedral_helix.h index 436895c5c3..172a8c3469 100644 --- a/src/EXTRA-MOLECULE/dihedral_helix.h +++ b/src/EXTRA-MOLECULE/dihedral_helix.h @@ -33,6 +33,7 @@ class DihedralHelix : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; protected: double *aphi, *bphi, *cphi; From 7ab9da0212f1f3601f7238fe2997f9313d8cc923 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 21 Jun 2023 10:59:13 +0300 Subject: [PATCH 18/60] Implementation of born_matrix in dihedral_helix.cpp --- src/EXTRA-MOLECULE/dihedral_helix.cpp | 106 ++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) diff --git a/src/EXTRA-MOLECULE/dihedral_helix.cpp b/src/EXTRA-MOLECULE/dihedral_helix.cpp index 059bef74a4..1d99de6ba9 100644 --- a/src/EXTRA-MOLECULE/dihedral_helix.cpp +++ b/src/EXTRA-MOLECULE/dihedral_helix.cpp @@ -41,6 +41,7 @@ using namespace MathConst; DihedralHelix::DihedralHelix(LAMMPS *lmp) : Dihedral(lmp) { writedata = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -324,3 +325,108 @@ void DihedralHelix::write_data(FILE *fp) for (int i = 1; i <= atom->ndihedraltypes; i++) fprintf(fp,"%d %g %g %g\n",i,aphi[i],bphi[i],cphi[i]); } + +/* ----------------------------------------------------------------------*/ + +void DihedralHelix::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double sb1,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; + double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; + double c2mag,sc1,sc2,s12,c; + double cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + + int **dihedrallist = neighbor->dihedrallist; + double **x = atom->x; + + int type = dihedrallist[nd][4]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + cx = vb1y*vb2z - vb1z*vb2y; + cy = vb1z*vb2x - vb1x*vb2z; + cz = vb1x*vb2y - vb1y*vb2x; + cmag = sqrt(cx*cx + cy*cy + cz*cz); + dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + phi = acos(c); + if (dx > 0.0) phi *= -1.0; + si = sin(phi); + if (fabs(si) < SMALLER) si = SMALLER; + siinv = 1.0/si; + + du = -aphi[type] + 3.0*bphi[type]*sin(3.0*phi)*siinv + + cphi[type]*sin(phi + MY_PI4)*siinv; + du2 = -(9.0*bphi[type]*cos(3.0*phi) + cphi[type]*cos(phi + MY_PI4))*siinv*siinv + + (3.0*bphi[type]*sin(3.0*phi) + cphi[type]*sin(phi + MY_PI4))*c*siinv*siinv*siinv; +} From ae96c9bd47e52d3c3066b40914b8e27ab6b042eb Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 22 Jun 2023 16:58:00 +0300 Subject: [PATCH 19/60] Define born_matrix() in dihedral_quadratic.h --- src/EXTRA-MOLECULE/dihedral_quadratic.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-MOLECULE/dihedral_quadratic.h b/src/EXTRA-MOLECULE/dihedral_quadratic.h index 90d8c3be6e..89f6fa3b25 100644 --- a/src/EXTRA-MOLECULE/dihedral_quadratic.h +++ b/src/EXTRA-MOLECULE/dihedral_quadratic.h @@ -33,6 +33,7 @@ class DihedralQuadratic : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; protected: double *k, *phi0; From 8e1711c8031699756ba3b5d84fa4582ad161a357 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 22 Jun 2023 17:00:10 +0300 Subject: [PATCH 20/60] Implement born_matrix in dihedral_quadratic.cpp --- src/EXTRA-MOLECULE/dihedral_quadratic.cpp | 110 ++++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/src/EXTRA-MOLECULE/dihedral_quadratic.cpp b/src/EXTRA-MOLECULE/dihedral_quadratic.cpp index cbe9e3e3a2..f576e6efdd 100644 --- a/src/EXTRA-MOLECULE/dihedral_quadratic.cpp +++ b/src/EXTRA-MOLECULE/dihedral_quadratic.cpp @@ -41,6 +41,7 @@ using namespace MathConst; DihedralQuadratic::DihedralQuadratic(LAMMPS *lmp) : Dihedral(lmp) { writedata = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -327,3 +328,112 @@ void DihedralQuadratic::write_data(FILE *fp) for (int i = 1; i <= atom->ndihedraltypes; i++) fprintf(fp,"%d %g %g \n",i,k[i],phi0[i]*180.0/MY_PI); } + +/* ----------------------------------------------------------------------*/ + +void DihedralQuadratic::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double sb1,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; + double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; + double c2mag,sc1,sc2,s12,c; + double s1,s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + + int **dihedrallist = neighbor->dihedrallist; + double **x = atom->x; + + int type = dihedrallist[nd][4]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + cx = vb1y*vb2z - vb1z*vb2y; + cy = vb1z*vb2x - vb1x*vb2z; + cz = vb1x*vb2y - vb1y*vb2x; + cmag = sqrt(cx*cx + cy*cy + cz*cz); + dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) + problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + phi = acos(c); + if (dx > 0.0) phi *= -1.0; + si = sin(phi); + if (fabs(si) < SMALLER) si = SMALLER; + siinv = 1.0/si; + + double dphi = phi-phi0[type]; + if (dphi > MY_PI) dphi -= 2*MY_PI; + else if (dphi < -MY_PI) dphi += 2*MY_PI; + + du = - 2.0 * k[type] * dphi * siinv; + du2 = 2.0 * k[type] * siinv * siinv * ( 1.0 - dphi * c * siinv) ; +} From 6e32d2932275271ced51152b3c849318c8956601 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 29 Jun 2023 08:53:50 +0300 Subject: [PATCH 21/60] clean up of originflag variable in compute_stress_mop_profile.cpp --- .../compute_stress_mop_profile.cpp | 38 ++++++++----------- 1 file changed, 16 insertions(+), 22 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 3a85869f3c..d37a941fde 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -38,7 +38,6 @@ using namespace LAMMPS_NS; enum { X, Y, Z }; -enum { LOWER, CENTER, UPPER, COORD }; enum { TOTAL, CONF, KIN, PAIR, BOND }; // clang-format off @@ -63,11 +62,15 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a // bin parameters - if (strcmp(arg[4],"lower") == 0) originflag = LOWER; - else if (strcmp(arg[4],"center") == 0) originflag = CENTER; - else if (strcmp(arg[4],"upper") == 0) originflag = UPPER; - else originflag = COORD; - if (originflag == COORD) origin = utils::numeric(FLERR,arg[4],false,lmp); + if (strcmp(arg[4],"lower") == 0) { + origin = domain->boxlo[dir]; + } else if (strcmp(arg[4],"center") == 0) { + origin = 0.5 * (domain->boxlo[dir] + domain->boxhi[dir]); + } else if (strcmp(arg[4],"upper") == 0) { + origin = domain->boxhi[dir]; + } else { + origin = utils::numeric(FLERR,arg[4],false,lmp); + } delta = utils::numeric(FLERR,arg[5],false,lmp); invdelta = 1.0/delta; @@ -620,23 +623,14 @@ void ComputeStressMopProfile::setup_bins() boxlo = domain->boxlo; boxhi = domain->boxhi; - if (originflag == LOWER) origin = boxlo[dir]; - else if (originflag == UPPER) origin = boxhi[dir]; - else if (originflag == CENTER) - origin = 0.5 * (boxlo[dir] + boxhi[dir]); + if ((origin > domain->boxhi[dir]) || (origin < domain->boxlo[dir])) + error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" ); - if (origin < boxlo[dir]) { - error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" ); - } else { - n = static_cast ((origin - boxlo[dir]) * invdelta); - lo = origin - n*delta; - } - if (origin < boxhi[dir]) { - n = static_cast ((boxhi[dir] - origin) * invdelta); - hi = origin + n*delta; - } else { - error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" ); - } + n = static_cast ((origin - boxlo[dir]) * invdelta); + lo = origin - n*delta; + + n = static_cast ((boxhi[dir] - origin) * invdelta); + hi = origin + n*delta; offset = lo; nbins = static_cast ((hi-lo) * invdelta + 1.5); From 484e7ad0e3c295514f063e1e7933f3683e93f383 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 29 Jun 2023 08:54:30 +0300 Subject: [PATCH 22/60] clean up of originflag from compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index 1b5b3a911f..5890505d71 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -46,7 +46,6 @@ class ComputeStressMopProfile : public Compute { int bondflag; - int originflag; double origin, delta, offset, invdelta; int nbins; double *coord, *coordp; From 005c15c07b3d4dd445ba276e8f3615852fdb7e71 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Jul 2023 13:26:47 +0300 Subject: [PATCH 23/60] compute kinetic contribution without assuming orthogonal geometry --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 98e3bf7043..18e6ec6115 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -422,6 +422,11 @@ void ComputeStressMop::compute_pairs() xi[1] = atom->x[i][1]; xi[2] = atom->x[i][2]; + // minimum image of xi with respect to the plane + xi[dir] -= pos; + domain->minimum_image(xi[0], xi[1], xi[2]); + xi[dir] += pos; + //velocities at t vi[0] = atom->v[i][0]; vi[1] = atom->v[i][1]; @@ -447,10 +452,8 @@ void ComputeStressMop::compute_pairs() // at each timestep, must check atoms going through the // image of the plane that is closest to the box - double pos_temp = pos+copysign(1.0,domain->prd_half[dir]-pos)*domain->prd[dir]; - if (fabs(xi[dir]-pos)= 0)) { // sgn = copysign(1.0,vi[dir]-vcm[dir]); sgn = copysign(1.0,vi[dir]); From 94fa2f51c9189874ab84661519f0047cd1f0c9db Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Jul 2023 13:38:56 +0300 Subject: [PATCH 24/60] compute kinetic contribution without assuming orthogonal geometry --- .../compute_stress_mop_profile.cpp | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index d37a941fde..c8793586ea 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -453,7 +453,22 @@ void ComputeStressMopProfile::compute_pairs() pos = coord[ibin]; pos1 = coordp[ibin]; - if (((xi[dir]-pos)*(xj[dir]-pos)*(xi[dir]-pos1)*(xj[dir]-pos1)<0)) { + // minimum image of xi with respect to the plane + xi[dir] -= pos; + domain->minimum_image(xi[0], xi[1], xi[2]); + xi[dir] += pos; + + // minimum image of xj with respect to xi + xj[0] -= xi[0]; + xj[1] -= xi[1]; + xj[2] -= xi[2]; + domain->minimum_image(xi[0], xi[1], xi[2]); + xj[0] += xi[0]; + xj[1] += xi[1]; + xj[2] += xi[2]; + + double tau = (xi[dir] - pos) / (xi[dir] - xj[dir]); + if ((tau <= 1) && (tau >= 0)) { sgn = copysign(1.0,vi[dir]); From 79ed2d9e8bb1bd5fefbf6fcd24eff5fd2d452381 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Jul 2023 16:35:25 +0300 Subject: [PATCH 25/60] Definition of compute_angle and related variables in compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index 5890505d71..c7214c055c 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -39,18 +39,20 @@ class ComputeStressMopProfile : public Compute { private: void compute_pairs(); void compute_bonds(); + void compute_angles(); void setup_bins(); int nvalues, dir; int *which; - int bondflag; + int bondflag, angleflag; double origin, delta, offset, invdelta; int nbins; double *coord, *coordp; double **values_local, **values_global; double **bond_local, **bond_global; + double **angle_local, **angle_global; double **local_contribution; double dt, nktv2p, ftm2v; From 9aa9bdd3ba4bd6f59cbe2723bedc46936cf77b62 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Jul 2023 16:45:53 +0300 Subject: [PATCH 26/60] Implementation of compute_angles in compute_stress_mop_profile.cpp and related adjustments to flags/memory allocations --- .../compute_stress_mop_profile.cpp | 227 +++++++++++++++++- 1 file changed, 221 insertions(+), 6 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index c8793586ea..dabbbc4781 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -13,7 +13,7 @@ /*------------------------------------------------------------------------ Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon) - Support for bonds added by : Evangelos Voyiatzis (NovaMechanics) + Support for bonds and angles added by : Evangelos Voyiatzis (NovaMechanics) --------------------------------------------------------------------------*/ #include "compute_stress_mop_profile.h" @@ -38,7 +38,7 @@ using namespace LAMMPS_NS; enum { X, Y, Z }; -enum { TOTAL, CONF, KIN, PAIR, BOND }; +enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE }; // clang-format off /* ---------------------------------------------------------------------- */ @@ -107,6 +107,11 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a which[nvalues] = BOND; nvalues++; } + } else if (strcmp(arg[iarg],"angle") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = ANGLE; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break; iarg++; @@ -131,6 +136,8 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a values_local = values_global = array = nullptr; bond_local = nullptr; bond_global = nullptr; + angle_local = nullptr; + angle_global = nullptr; local_contribution = nullptr; // bin setup @@ -159,6 +166,8 @@ ComputeStressMopProfile::~ComputeStressMopProfile() memory->destroy(values_global); memory->destroy(bond_local); memory->destroy(bond_global); + memory->destroy(angle_local); + memory->destroy(angle_global); memory->destroy(local_contribution); memory->destroy(array); } @@ -206,9 +215,14 @@ void ComputeStressMopProfile::init() if (force->bond) bondflag = 1; - if (force->angle) - if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) - error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); + if (force->angle) { + if (force->angle->born_matrix_enable == 0) { + if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) + error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); + } else { + angleflag = 1; + } + } if (force->dihedral) if ((strcmp(force->dihedral_style, "zero") != 0) && (strcmp(force->dihedral_style, "none") != 0)) error->all(FLERR,"compute stress/mop/profile does not account for dihedral potentials"); @@ -260,13 +274,27 @@ void ComputeStressMopProfile::compute_array() // sum bond contribution over all procs MPI_Allreduce(&bond_local[0][0],&bond_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + if (angleflag) { + //Compute angle contribution on separate procs + compute_angles(); + } else { + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + angle_local[m][i] = 0.0; + } + } + } + + // sum angle contribution over all procs + MPI_Allreduce(&angle_local[0][0],&angle_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + for (int ibin=0; ibinx; + tagint *tag = atom->tag; + int *num_angle = atom->num_angle; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their angles + Angle *angle = force->angle; + + double duang, du2ang; + double dx[3] = {0.0, 0.0, 0.0}; + double dx_left[3] = {0.0, 0.0, 0.0}; + double dx_right[3] = {0.0, 0.0, 0.0}; + double x_angle_left[3] = {0.0, 0.0, 0.0}; + double x_angle_middle[3] = {0.0, 0.0, 0.0}; + double x_angle_right[3] = {0.0, 0.0, 0.0}; + double dcos_theta[3] = {0.0, 0.0, 0.0}; + + // initialization + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + angle_local[m][i] = 0.0; + } + local_contribution[m][0] = 0.0; + local_contribution[m][1] = 0.0; + local_contribution[m][2] = 0.0; + } + + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == 1) + na = num_angle[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + na = onemols[imol]->num_angle[iatom]; + } + + for (int i = 0; i < na; i++) { + if (molecular == 1) { + if (tag[atom2] != angle_atom2[atom2][i]) continue; + atype = angle_type[atom2][i]; + atom1 = atom->map(angle_atom1[atom2][i]); + atom3 = atom->map(angle_atom3[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; + atype = onemols[imol]->angle_type[atom2][i]; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atype <= 0) continue; + + for (int ibin = 0; ibinminimum_image(dx[0], dx[1], dx[2]); + x_angle_left[0] = dx[0]; + x_angle_left[1] = dx[1]; + x_angle_left[2] = dx[2]; + x_angle_left[dir] += pos; + + // minimum image of atom2 with respect to atom1 + dx_left[0] = x[atom2][0] - x_angle_left[0]; + dx_left[1] = x[atom2][1] - x_angle_left[1]; + dx_left[2] = x[atom2][2] - x_angle_left[2]; + domain->minimum_image(dx_left[0], dx_left[1], dx_left[2]); + x_angle_middle[0] = x_angle_left[0] + dx_left[0]; + x_angle_middle[1] = x_angle_left[1] + dx_left[1]; + x_angle_middle[2] = x_angle_left[2] + dx_left[2]; + + // minimum image of atom3 with respect to atom2 + dx_right[0] = x[atom3][0] - x_angle_middle[0]; + dx_right[1] = x[atom3][1] - x_angle_middle[1]; + dx_right[2] = x[atom3][2] - x_angle_middle[2]; + domain->minimum_image(dx_right[0], dx_right[1], dx_right[2]); + x_angle_right[0] = x_angle_middle[0] + dx_right[0]; + x_angle_right[1] = x_angle_middle[1] + dx_right[1]; + x_angle_right[2] = x_angle_middle[2] + dx_right[2]; + + // check if any bond vector crosses the plane of interest + double tau_right = (x_angle_right[dir] - pos) / (x_angle_right[dir] - x_angle_middle[dir]); + double tau_left = (x_angle_middle[dir] - pos) / (x_angle_middle[dir] - x_angle_left[dir]); + bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + + // no bonds crossing the plane + if (!right_cross && !left_cross) continue; + + // compute the cos(theta) of the angle + r1 = sqrt(dx_left[0]*dx_left[0] + dx_left[1]*dx_left[1] + dx_left[2]*dx_left[2]); + r2 = sqrt(dx_right[0]*dx_right[0] + dx_right[1]*dx_right[1] + dx_right[2]*dx_right[2]); + cos_theta = -(dx_right[0]*dx_left[0] + dx_right[1]*dx_left[1] + dx_right[2]*dx_left[2])/(r1*r2); + + if (cos_theta > 1.0) cos_theta = 1.0; + if (cos_theta < -1.0) cos_theta = -1.0; + + // The method returns derivative with regards to cos(theta) + angle->born_matrix(atype, atom1, atom2, atom3, duang, du2ang); + // only right bond crossing the plane + if (right_cross && !left_cross) + { + double sgn = copysign(1.0, x_angle_right[dir] - pos); + dcos_theta[0] = sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2; + dcos_theta[1] = sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2; + dcos_theta[2] = sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2; + } + + // only left bond crossing the plane + if (!right_cross && left_cross) + { + double sgn = copysign(1.0, x_angle_left[dir] - pos); + dcos_theta[0] = -sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1; + dcos_theta[1] = -sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1; + dcos_theta[2] = -sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1; + } + + // both bonds crossing the plane + if (right_cross && left_cross) + { + // due to right bond + double sgn = copysign(1.0, x_angle_middle[dir] - pos); + dcos_theta[0] = -sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2; + dcos_theta[1] = -sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2; + dcos_theta[2] = -sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2; + + // due to left bond + dcos_theta[0] += sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1; + dcos_theta[1] += sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1; + dcos_theta[2] += sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1; + } + + // final contribution of the given angle term + local_contribution[ibin][0] += duang*dcos_theta[0]/area*nktv2p; + local_contribution[ibin][1] += duang*dcos_theta[1]/area*nktv2p; + local_contribution[ibin][2] += duang*dcos_theta[2]/area*nktv2p; + } + } + } + + // loop over the keywords and if necessary add the angle contribution + int m = 0; + while (m < nvalues) { + if (which[m] == CONF || which[m] == TOTAL || which[m] == ANGLE) { + for (int ibin = 0; ibin < nbins; ibin++) { + angle_local[ibin][m] = local_contribution[ibin][0]; + angle_local[ibin][m+1] = local_contribution[ibin][1]; + angle_local[ibin][m+2] = local_contribution[ibin][2]; + } + } + m += 3; + } +} + /* ---------------------------------------------------------------------- setup 1d bins and their extent and coordinates called at init() @@ -657,6 +870,8 @@ void ComputeStressMopProfile::setup_bins() memory->create(values_global,nbins,nvalues,"stress/mop/profile:values_global"); memory->create(bond_local,nbins,nvalues,"stress/mop/profile:bond_local"); memory->create(bond_global,nbins,nvalues,"stress/mop/profile:bond_global"); + memory->create(angle_local,nbins,nvalues,"stress/mop/profile:angle_local"); + memory->create(angle_global,nbins,nvalues,"stress/mop/profile:angle_global"); memory->create(local_contribution,nbins,3,"stress/mop/profile:local_contribution"); // set bin coordinates From 78f4e4f1a15061991a7a0073566df2ed43712715 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Jul 2023 16:51:41 +0300 Subject: [PATCH 27/60] Update compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 21c2963545..70986862fe 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -74,9 +74,7 @@ Between one and six keywords can be used to indicate which contributions to the stress must be computed: total stress (total), kinetic stress (kin), configurational stress (conf), stress due to bond stretching (bond), stress due to angle bending (angle) and/or due to pairwise -non-bonded interactions (pair). The angle keyword is currently -available only for the *stress/mop* command and **not** the -*stress/mop/profile* command. +non-bonded interactions (pair). NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. @@ -136,7 +134,7 @@ requires the class method ``Pair::single()`` to be implemented, which is not possible for manybody potentials. In particular, compute *stress/mop/profile* does not work with more than two-body pair interactions, long range (kspace) interactions and -angle/dihedral/improper intramolecular interactions. Similarly, compute +dihedral/improper intramolecular interactions. Similarly, compute *stress/mop* does not work with more than two-body pair interactions, long range (kspace) interactions and dihedral/improper intramolecular interactions but works with all bond interactions with the class method From 9fde61fc4e358e44854cf73f8d5cb230b3bf7e7c Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Jul 2023 16:59:02 +0300 Subject: [PATCH 28/60] Update compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index dabbbc4781..2dee2bf60f 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -18,6 +18,7 @@ #include "compute_stress_mop_profile.h" +#include "angle.h" #include "atom.h" #include "atom_vec.h" #include "bond.h" From b81df7c21b6be9739f234d4dc9c73272b38cde9c Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:12:09 +0300 Subject: [PATCH 29/60] Include methods and variables for dihedral contribution to compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index c7214c055c..b9b97617c0 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -40,12 +40,13 @@ class ComputeStressMopProfile : public Compute { void compute_pairs(); void compute_bonds(); void compute_angles(); + void compute_dihedrals(); void setup_bins(); int nvalues, dir; int *which; - int bondflag, angleflag; + int bondflag, angleflag, dihedralflag; double origin, delta, offset, invdelta; int nbins; @@ -53,6 +54,7 @@ class ComputeStressMopProfile : public Compute { double **values_local, **values_global; double **bond_local, **bond_global; double **angle_local, **angle_global; + double **dihedral_local, **dihedral_global; double **local_contribution; double dt, nktv2p, ftm2v; From e819b38a1831a56839cc13be2877099244ca4afe Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:15:12 +0300 Subject: [PATCH 30/60] undo dihedral changes in compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index b9b97617c0..c7214c055c 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -40,13 +40,12 @@ class ComputeStressMopProfile : public Compute { void compute_pairs(); void compute_bonds(); void compute_angles(); - void compute_dihedrals(); void setup_bins(); int nvalues, dir; int *which; - int bondflag, angleflag, dihedralflag; + int bondflag, angleflag; double origin, delta, offset, invdelta; int nbins; @@ -54,7 +53,6 @@ class ComputeStressMopProfile : public Compute { double **values_local, **values_global; double **bond_local, **bond_global; double **angle_local, **angle_global; - double **dihedral_local, **dihedral_global; double **local_contribution; double dt, nktv2p, ftm2v; From 381d8de017c37794682e77a3825f4e1b9ddd4458 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:22:09 +0300 Subject: [PATCH 31/60] initialization of angleflag in constructor compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 2a82cc0f8b..3dfc9326ac 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -49,6 +49,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a if (narg < 7) utils::missing_cmd_args(FLERR, "compute stress/mop/profile", error); bondflag = 0; + angleflag = 0; // set compute mode and direction of plane(s) for pressure calculation From ec458e2861b5b4950bfcd9a9dfb67bf1c636c319 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:32:26 +0300 Subject: [PATCH 32/60] method and member variables definition for dihedrals in compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index c7214c055c..b9b97617c0 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -40,12 +40,13 @@ class ComputeStressMopProfile : public Compute { void compute_pairs(); void compute_bonds(); void compute_angles(); + void compute_dihedrals(); void setup_bins(); int nvalues, dir; int *which; - int bondflag, angleflag; + int bondflag, angleflag, dihedralflag; double origin, delta, offset, invdelta; int nbins; @@ -53,6 +54,7 @@ class ComputeStressMopProfile : public Compute { double **values_local, **values_global; double **bond_local, **bond_global; double **angle_local, **angle_global; + double **dihedral_local, **dihedral_global; double **local_contribution; double dt, nktv2p, ftm2v; From d40fb4a337c67c2f0c3bc79bf29b87353cd80ceb Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:45:46 +0300 Subject: [PATCH 33/60] method implementation for dihedral contribution to compute_stress_mop_profile.cpp --- .../compute_stress_mop_profile.cpp | 361 +++++++++++++++++- 1 file changed, 354 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 3dfc9326ac..27e420e71d 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -13,7 +13,7 @@ /*------------------------------------------------------------------------ Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon) - Support for bonds and angles added by : Evangelos Voyiatzis (NovaMechanics) + Support for bonds, angles and dihedrals added by : Evangelos Voyiatzis (NovaMechanics) --------------------------------------------------------------------------*/ #include "compute_stress_mop_profile.h" @@ -23,6 +23,7 @@ #include "atom_vec.h" #include "bond.h" #include "comm.h" +#include "dihedral.h" #include "domain.h" #include "error.h" #include "force.h" @@ -38,8 +39,10 @@ using namespace LAMMPS_NS; +#define SMALL 0.001 + enum { X, Y, Z }; -enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE }; +enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; /* ---------------------------------------------------------------------- */ @@ -50,6 +53,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a bondflag = 0; angleflag = 0; + dihedralflag = 0; // set compute mode and direction of plane(s) for pressure calculation @@ -114,6 +118,11 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a which[nvalues] = ANGLE; nvalues++; } + } else if (strcmp(arg[iarg],"dihedral") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = DIHEDRAL; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break; @@ -141,6 +150,8 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a bond_global = nullptr; angle_local = nullptr; angle_global = nullptr; + dihedral_local = nullptr; + dihedral_global = nullptr; local_contribution = nullptr; // bin setup @@ -171,6 +182,8 @@ ComputeStressMopProfile::~ComputeStressMopProfile() memory->destroy(bond_global); memory->destroy(angle_local); memory->destroy(angle_global); + memory->destroy(dihedral_local); + memory->destroy(dihedral_global); memory->destroy(local_contribution); memory->destroy(array); } @@ -227,10 +240,16 @@ void ComputeStressMopProfile::init() } } - if (force->dihedral) - if ((strcmp(force->dihedral_style, "zero") != 0) && - (strcmp(force->dihedral_style, "none") != 0)) - error->all(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); + if (force->dihedral) { + if (force->dihedral->born_matrix_enable == 0) { + if ((strcmp(force->dihedral_style, "zero") != 0) && + (strcmp(force->dihedral_style, "none") != 0)) + error->all(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); + } else { + dihedralflag = 1; + } + } + if (force->improper) if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) @@ -294,6 +313,17 @@ void ComputeStressMopProfile::compute_array() // sum angle contribution over all procs MPI_Allreduce(&angle_local[0][0],&angle_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + + if (dihedralflag) { + //Compute dihedral contribution on separate procs + compute_dihedrals(); + } else { + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + dihedral_local[m][i] = 0.0; + } + } + } for (int ibin = 0; ibin < nbins; ibin++) { array[ibin][0] = coord[ibin]; @@ -301,7 +331,7 @@ void ComputeStressMopProfile::compute_array() int mo = 1; int m = 0; while (m < nvalues) { - array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m]; + array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] + dihedral_global[ibin][m]; m++; } } @@ -835,6 +865,321 @@ void ComputeStressMopProfile::compute_angles() } } +/*------------------------------------------------------------------------ + compute dihedral contribution to pressure of local proc + -------------------------------------------------------------------------*/ + +void ComputeStressMopProfile::compute_dihedrals() +{ + int i, nd, atom1, atom2, atom3, atom4, imol, iatom; + tagint tagprev; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z; + double vb2xm, vb2ym, vb2zm; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, c2mag, ctmp, r12c1, c1mag, r12c2; + double s1, s2, s12, sc1, sc2, a11, a22, a33, a12, a13, a23; + double df[3], f1[3], f2[3], f3[3], f4[3]; + double c, sx2, sy2, sz2, sin2; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their dihedrals + + Dihedral *dihedral = force->dihedral; + + double dudih, du2dih; + + double diffx[3] = {0.0, 0.0, 0.0}; + double x_atom_1[3] = {0.0, 0.0, 0.0}; + double x_atom_2[3] = {0.0, 0.0, 0.0}; + double x_atom_3[3] = {0.0, 0.0, 0.0}; + double x_atom_4[3] = {0.0, 0.0, 0.0}; + + // initialization + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + dihedral_local[m][i] = 0.0; + } + local_contribution[m][0] = 0.0; + local_contribution[m][1] = 0.0; + local_contribution[m][2] = 0.0; + } + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == Atom::MOLECULAR) + nd = num_dihedral[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + nd = onemols[imol]->num_dihedral[iatom]; + } + + for (i = 0; i < nd; i++) { + if (molecular == 1) { + if (tag[atom2] != dihedral_atom2[atom2][i]) continue; + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; + + for (int ibin = 0; ibinminimum_image(x_atom_1[0], x_atom_1[1], x_atom_1[2]); + x_atom_1[dir] += pos; + + // minimum image of atom2 with respect to atom1 + diffx[0] = x[atom2][0] - x_atom_1[0]; + diffx[1] = x[atom2][1] - x_atom_1[1]; + diffx[2] = x[atom2][2] - x_atom_1[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_2[0] = x_atom_1[0] + diffx[0]; + x_atom_2[1] = x_atom_1[1] + diffx[1]; + x_atom_2[2] = x_atom_1[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom3][0] - x_atom_2[0]; + diffx[1] = x[atom3][1] - x_atom_2[1]; + diffx[2] = x[atom3][2] - x_atom_2[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_3[0] = x_atom_2[0] + diffx[0]; + x_atom_3[1] = x_atom_2[1] + diffx[1]; + x_atom_3[2] = x_atom_2[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom4][0] - x_atom_3[0]; + diffx[1] = x[atom4][1] - x_atom_3[1]; + diffx[2] = x[atom4][2] - x_atom_3[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_4[0] = x_atom_3[0] + diffx[0]; + x_atom_4[1] = x_atom_3[1] + diffx[1]; + x_atom_4[2] = x_atom_3[2] + diffx[2]; + + // check if any bond vector crosses the plane of interest + double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]); + double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]); + double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]); + bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1)); + bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + + // no bonds crossing the plane + if (!right_cross && !middle_cross && !left_cross) continue; + + dihedral->born_matrix(i, atom1, atom2, atom3, atom4, dudih, du2dih); + + // first bond + vb1x = x_atom_1[0] - x_atom_2[0]; + vb1y = x_atom_1[1] - x_atom_2[1]; + vb1z = x_atom_1[2] - x_atom_2[2]; + + // second bond + vb2x = x_atom_3[0] - x_atom_2[0]; + vb2y = x_atom_3[1] - x_atom_2[1]; + vb2z = x_atom_3[2] - x_atom_2[2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // third bond + vb3x = x_atom_4[0] - x_atom_3[0]; + vb3y = x_atom_4[1] - x_atom_3[1]; + vb3z = x_atom_4[2] - x_atom_3[2]; + + // c0 calculation + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + // 1st and 2nd angle + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + // error check + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // forces on each particle + double a = dudih; + c = c * a; + s12 = s12 * a; + a11 = c*sb1*s1; + a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); + a33 = c*sb3*s2; + a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); + a13 = -rb1*rb3*s12; + a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + + sx2 = a12*vb1x + a22*vb2x + a23*vb3x; + sy2 = a12*vb1y + a22*vb2y + a23*vb3y; + sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + + f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; + f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; + f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + + f2[0] = -sx2 - f1[0]; + f2[1] = -sy2 - f1[1]; + f2[2] = -sz2 - f1[2]; + + f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; + f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; + f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // only right bond crossing the plane + if (right_cross && !middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * f1[0]; + df[1] = sgn * f1[1]; + df[2] = sgn * f1[2]; + } + + // only middle bond crossing the plane + if (!right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * (f2[0] + f1[0]); + df[1] = sgn * (f2[1] + f1[1]); + df[2] = sgn * (f2[2] + f1[2]); + } + + // only left bond crossing the plane + if (!right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_4[dir] - pos); + df[0] = sgn * f4[0]; + df[1] = sgn * f4[1]; + df[2] = sgn * f4[2]; + } + + // only right & middle bonds crossing the plane + if (right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * f2[0]; + df[1] = sgn * f2[1]; + df[2] = sgn * f2[2]; + } + + // only right & left bonds crossing the plane + if (right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f4[0]); + df[1] = sgn * (f1[1] + f4[1]); + df[2] = sgn * (f1[2] + f4[2]); + } + + // only middle & left bonds crossing the plane + if (!right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_3[dir] - pos); + df[0] = sgn * f3[0]; + df[1] = sgn * f3[1]; + df[2] = sgn * f3[2]; + } + + // all three bonds crossing the plane + if (right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f3[0]); + df[1] = sgn * (f1[1] + f3[1]); + df[2] = sgn * (f1[2] + f3[2]); + } + + local_contribution[ibin][0] += df[0]/area*nktv2p; + local_contribution[ibin][1] += df[1]/area*nktv2p; + local_contribution[ibin][2] += df[2]/area*nktv2p; + } + } + } + + // loop over the keywords and if necessary add the dihedral contribution + int m = 0; + while (m < nvalues) { + if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) { + for (int ibin = 0; ibin < nbins; ibin++) { + dihedral_local[ibin][m] = local_contribution[ibin][0]; + dihedral_local[ibin][m+1] = local_contribution[ibin][1]; + dihedral_local[ibin][m+2] = local_contribution[ibin][2]; + } + } + m += 3; + } + +} + /* ---------------------------------------------------------------------- setup 1d bins and their extent and coordinates called at init() @@ -870,6 +1215,8 @@ void ComputeStressMopProfile::setup_bins() memory->create(bond_global, nbins, nvalues, "stress/mop/profile:bond_global"); memory->create(angle_local, nbins, nvalues, "stress/mop/profile:angle_local"); memory->create(angle_global, nbins, nvalues, "stress/mop/profile:angle_global"); + memory->create(dihedral_local,nbins,nvalues,"stress/mop/profile:dihedral_local"); + memory->create(dihedral_global,nbins,nvalues,"stress/mop/profile:dihedral_global"); memory->create(local_contribution, nbins, 3, "stress/mop/profile:local_contribution"); // set bin coordinates From b86d1f655381c076a6b378bc8ea8f46aedf5c6de Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:48:45 +0300 Subject: [PATCH 34/60] Update compute_stress_mop.rst for dihedral interactions --- doc/src/compute_stress_mop.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 70986862fe..3cd7a67c9c 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -18,7 +18,7 @@ Syntax * style = *stress/mop* or *stress/mop/profile* * dir = *x* or *y* or *z* is the direction normal to the plane * args = argument specific to the compute style -* keywords = *kin* or *conf* or *total* or *pair* or *bond* or *angle* (one or more can be specified) +* keywords = *kin* or *conf* or *total* or *pair* or *bond* or *angle* or *dihedral* (one or more can be specified) .. parsed-literal:: @@ -68,13 +68,13 @@ Verlet algorithm. .. versionadded:: 15Jun2023 - contributions from bond and angle potentials + contributions from bond, angle and dihedral potentials -Between one and six keywords can be used to indicate which contributions +Between one and seven keywords can be used to indicate which contributions to the stress must be computed: total stress (total), kinetic stress (kin), configurational stress (conf), stress due to bond stretching -(bond), stress due to angle bending (angle) and/or due to pairwise -non-bonded interactions (pair). +(bond), stress due to angle bending (angle), stress due to dihedral terms (dihedral) +and/or due to pairwise non-bonded interactions (pair). NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. @@ -134,7 +134,7 @@ requires the class method ``Pair::single()`` to be implemented, which is not possible for manybody potentials. In particular, compute *stress/mop/profile* does not work with more than two-body pair interactions, long range (kspace) interactions and -dihedral/improper intramolecular interactions. Similarly, compute +improper intramolecular interactions. Similarly, compute *stress/mop* does not work with more than two-body pair interactions, long range (kspace) interactions and dihedral/improper intramolecular interactions but works with all bond interactions with the class method From 1591b21617a196607306955e7d108f38b330ed1d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 26 Sep 2023 16:53:16 +0300 Subject: [PATCH 35/60] remove whitespaces from compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 27e420e71d..223e4da66a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -324,7 +324,7 @@ void ComputeStressMopProfile::compute_array() } } } - + for (int ibin = 0; ibin < nbins; ibin++) { array[ibin][0] = coord[ibin]; @@ -540,7 +540,7 @@ void ComputeStressMopProfile::compute_pairs() vcross[0] = vi[0] - fi[0] * iterm; vcross[1] = vi[1] - fi[1] * iterm; vcross[2] = vi[2] - fi[2] * iterm; - + values_local[ibin][m] += imass * vcross[0] * sgn / dt / area * nktv2p / ftm2v; values_local[ibin][m + 1] += imass * vcross[1] * sgn / dt / area * nktv2p / ftm2v; values_local[ibin][m + 2] += imass * vcross[2] * sgn / dt / area * nktv2p / ftm2v; From 3445330cf1679addfdce5ab3e678b54669dc933d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 28 Sep 2023 16:13:13 +0300 Subject: [PATCH 36/60] remove whitespace from compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index a1077660e5..3f1ae008ea 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -433,7 +433,7 @@ void ComputeStressMop::compute_pairs() xi[dir] -= pos; domain->minimum_image(xi[0], xi[1], xi[2]); xi[dir] += pos; - + //velocities at t vi[0] = atom->v[i][0]; From ac435319fdf8b127ac33cfdf24ea6b8bfcc4649b Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 28 Sep 2023 16:15:03 +0300 Subject: [PATCH 37/60] Definition of compute_dihedral and related variables in compute_stress_mop.h --- src/EXTRA-COMPUTE/compute_stress_mop.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index 86140dc278..0a0ea8b55a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -40,15 +40,17 @@ class ComputeStressMop : public Compute { void compute_pairs(); void compute_bonds(); void compute_angles(); + void compute_dihedrals(); int nvalues, dir; int *which; - int bondflag, angleflag; + int bondflag, angleflag, dihedralflag; double *values_local, *values_global; double *bond_local, *bond_global; double *angle_local, *angle_global; + double *dihedral_local, *dihedral_global; double pos, pos1, dt, nktv2p, ftm2v; double area; class NeighList *list; From ca449f1ea859d47e665ae98e5a58cf8aa9781e60 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 28 Sep 2023 16:25:52 +0300 Subject: [PATCH 38/60] Prepare for inclusion of dihedral contribution in compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 47 +++++++++++++++++++++--- 1 file changed, 42 insertions(+), 5 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 3f1ae008ea..d94ba61d71 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -23,6 +23,7 @@ #include "atom_vec.h" #include "bond.h" #include "comm.h" +#include "dihedral.h" #include "domain.h" #include "error.h" #include "force.h" @@ -38,8 +39,10 @@ using namespace LAMMPS_NS; +#define SMALL 0.001 + enum { X, Y, Z }; -enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE }; +enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; /* ---------------------------------------------------------------------- */ @@ -49,6 +52,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( bondflag = 0; angleflag = 0; + dihedralflag = 0; // set compute mode and direction of plane(s) for pressure calculation @@ -129,6 +133,11 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( which[nvalues] = ANGLE; nvalues++; } + } else if (strcmp(arg[iarg],"dihedral") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = DIHEDRAL; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop command"); //break; @@ -152,6 +161,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( bond_global = nullptr; angle_local = nullptr; angle_global = nullptr; + dihedral_local = nullptr; + dihedral_global = nullptr; // this fix produces a global vector @@ -162,6 +173,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( memory->create(bond_global, nvalues, "stress/mop:bond_global"); memory->create(angle_local, nvalues, "stress/mop:angle_local"); memory->create(angle_global, nvalues, "stress/mop:angle_global"); + memory->create(dihedral_local,nvalues,"stress/mop:dihedral_local"); + memory->create(dihedral_global,nvalues,"stress/mop:dihedral_global"); size_vector = nvalues; vector_flag = 1; @@ -180,6 +193,8 @@ ComputeStressMop::~ComputeStressMop() memory->destroy(bond_global); memory->destroy(angle_local); memory->destroy(angle_global); + memory->destroy(dihedral_local); + memory->destroy(dihedral_global); memory->destroy(vector); } @@ -233,9 +248,13 @@ void ComputeStressMop::init() } } if (force->dihedral) { - if ((strcmp(force->dihedral_style, "zero") != 0) && - (strcmp(force->dihedral_style, "none") != 0)) - error->all(FLERR, "compute stress/mop does not account for dihedral potentials"); + if (force->dihedral->born_matrix_enable == 0) { + if ((strcmp(force->dihedral_style, "zero") != 0) && + (strcmp(force->dihedral_style, "none") != 0)) + error->all(FLERR, "compute stress/mop does not account for dihedral potentials"); + } else { + dihedralflag = 1; + } } if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && @@ -297,8 +316,18 @@ void ComputeStressMop::compute_vector() MPI_Allreduce(angle_local, angle_global, nvalues, MPI_DOUBLE, MPI_SUM, world); + if (dihedralflag) { + //Compute dihedral contribution on separate procs + compute_dihedrals(); + } else { + for (int i=0; i Date: Thu, 28 Sep 2023 16:29:15 +0300 Subject: [PATCH 39/60] remove whitespace from compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index d94ba61d71..331f4e4841 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -325,7 +325,7 @@ void ComputeStressMop::compute_vector() // sum dihedral contribution over all procs MPI_Allreduce(dihedral_local,dihedral_global,nvalues,MPI_DOUBLE,MPI_SUM,world); - + for (int m = 0; m < nvalues; m++) { vector[m] = values_global[m] + bond_global[m] + angle_global[m] + dihedral_global[m]; } From bbd6b2846f83cef48ed3960fe2496c64ee37f99e Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 28 Sep 2023 16:59:03 +0300 Subject: [PATCH 40/60] implementation of compute_dihedral() in compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 297 +++++++++++++++++++++++ 1 file changed, 297 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 331f4e4841..6c35b4ba07 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -825,4 +825,301 @@ void ComputeStressMop::compute_angles() void ComputeStressMop::compute_dihedrals() { + int i, nd, atom1, atom2, atom3, atom4, imol, iatom; + tagint tagprev; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z; + double vb2xm, vb2ym, vb2zm; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, c2mag, ctmp, r12c1, c1mag, r12c2; + double s1, s2, s12, sc1, sc2, a11, a22, a33, a12, a13, a23; + double df[3], f1[3], f2[3], f3[3], f4[3]; + double c, sx2, sy2, sz2, sin2; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their dihedrals + + Dihedral *dihedral = force->dihedral; + + double dudih, du2dih; + + double diffx[3] = {0.0, 0.0, 0.0}; + double x_atom_1[3] = {0.0, 0.0, 0.0}; + double x_atom_2[3] = {0.0, 0.0, 0.0}; + double x_atom_3[3] = {0.0, 0.0, 0.0}; + double x_atom_4[3] = {0.0, 0.0, 0.0}; + + // initialization + for (int i = 0; i < nvalues; i++) { + dihedral_local[i] = 0.0; + } + double local_contribution[3] = {0.0, 0.0, 0.0}; + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == Atom::MOLECULAR) + nd = num_dihedral[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + nd = onemols[imol]->num_dihedral[iatom]; + } + + for (i = 0; i < nd; i++) { + if (molecular == 1) { + if (tag[atom2] != dihedral_atom2[atom2][i]) continue; + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; + + // minimum image of atom1 with respect to the plane of interest + x_atom_1[0] = x[atom1][0]; + x_atom_1[1] = x[atom1][1]; + x_atom_1[2] = x[atom1][2]; + x_atom_1[dir] -= pos; + domain->minimum_image(x_atom_1[0], x_atom_1[1], x_atom_1[2]); + x_atom_1[dir] += pos; + + // minimum image of atom2 with respect to atom1 + diffx[0] = x[atom2][0] - x_atom_1[0]; + diffx[1] = x[atom2][1] - x_atom_1[1]; + diffx[2] = x[atom2][2] - x_atom_1[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_2[0] = x_atom_1[0] + diffx[0]; + x_atom_2[1] = x_atom_1[1] + diffx[1]; + x_atom_2[2] = x_atom_1[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom3][0] - x_atom_2[0]; + diffx[1] = x[atom3][1] - x_atom_2[1]; + diffx[2] = x[atom3][2] - x_atom_2[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_3[0] = x_atom_2[0] + diffx[0]; + x_atom_3[1] = x_atom_2[1] + diffx[1]; + x_atom_3[2] = x_atom_2[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom4][0] - x_atom_3[0]; + diffx[1] = x[atom4][1] - x_atom_3[1]; + diffx[2] = x[atom4][2] - x_atom_3[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_4[0] = x_atom_3[0] + diffx[0]; + x_atom_4[1] = x_atom_3[1] + diffx[1]; + x_atom_4[2] = x_atom_3[2] + diffx[2]; + + // check if any bond vector crosses the plane of interest + double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]); + double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]); + double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]); + bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1)); + bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + + // no bonds crossing the plane + if (!right_cross && !middle_cross && !left_cross) continue; + + dihedral->born_matrix(i, atom1, atom2, atom3, atom4, dudih, du2dih); + + // first bond + vb1x = x_atom_1[0] - x_atom_2[0]; + vb1y = x_atom_1[1] - x_atom_2[1]; + vb1z = x_atom_1[2] - x_atom_2[2]; + + // second bond + vb2x = x_atom_3[0] - x_atom_2[0]; + vb2y = x_atom_3[1] - x_atom_2[1]; + vb2z = x_atom_3[2] - x_atom_2[2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // third bond + vb3x = x_atom_4[0] - x_atom_3[0]; + vb3y = x_atom_4[1] - x_atom_3[1]; + vb3z = x_atom_4[2] - x_atom_3[2]; + + // c0 calculation + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + // 1st and 2nd angle + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + // error check + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // forces on each particle + double a = dudih; + c = c * a; + s12 = s12 * a; + a11 = c*sb1*s1; + a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); + a33 = c*sb3*s2; + a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); + a13 = -rb1*rb3*s12; + a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + + sx2 = a12*vb1x + a22*vb2x + a23*vb3x; + sy2 = a12*vb1y + a22*vb2y + a23*vb3y; + sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + + f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; + f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; + f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + + f2[0] = -sx2 - f1[0]; + f2[1] = -sy2 - f1[1]; + f2[2] = -sz2 - f1[2]; + + f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; + f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; + f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // only right bond crossing the plane + if (right_cross && !middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * f1[0]; + df[1] = sgn * f1[1]; + df[2] = sgn * f1[2]; + } + + // only middle bond crossing the plane + if (!right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * (f2[0] + f1[0]); + df[1] = sgn * (f2[1] + f1[1]); + df[2] = sgn * (f2[2] + f1[2]); + } + + // only left bond crossing the plane + if (!right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_4[dir] - pos); + df[0] = sgn * f4[0]; + df[1] = sgn * f4[1]; + df[2] = sgn * f4[2]; + } + + // only right & middle bonds crossing the plane + if (right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * f2[0]; + df[1] = sgn * f2[1]; + df[2] = sgn * f2[2]; + } + + // only right & left bonds crossing the plane + if (right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f4[0]); + df[1] = sgn * (f1[1] + f4[1]); + df[2] = sgn * (f1[2] + f4[2]); + } + + // only middle & left bonds crossing the plane + if (!right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_3[dir] - pos); + df[0] = sgn * f3[0]; + df[1] = sgn * f3[1]; + df[2] = sgn * f3[2]; + } + + // all three bonds crossing the plane + if (right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f3[0]); + df[1] = sgn * (f1[1] + f3[1]); + df[2] = sgn * (f1[2] + f3[2]); + } + + local_contribution[0] += df[0]/area*nktv2p; + local_contribution[1] += df[1]/area*nktv2p; + local_contribution[2] += df[2]/area*nktv2p; + } + } + + // loop over the keywords and if necessary add the dihedral contribution + int m = 0; + while (m < nvalues) { + if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) { + dihedral_local[m] = local_contribution[0]; + dihedral_local[m+1] = local_contribution[1]; + dihedral_local[m+2] = local_contribution[2]; + } + m += 3; + } + } From dc84ab5e5fc826e8095eb7d2c2ec38628dc4c7f7 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 28 Sep 2023 17:06:11 +0300 Subject: [PATCH 41/60] Update compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 3cd7a67c9c..74d4c618e7 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -132,14 +132,9 @@ size does not change in time, and axis-aligned planes. The method only works with two-body pair interactions, because it requires the class method ``Pair::single()`` to be implemented, which is not possible for manybody potentials. In particular, compute -*stress/mop/profile* does not work with more than two-body pair +*stress/mop/profile* and *stress/mop* do not work with more than two-body pair interactions, long range (kspace) interactions and -improper intramolecular interactions. Similarly, compute -*stress/mop* does not work with more than two-body pair interactions, -long range (kspace) interactions and dihedral/improper intramolecular -interactions but works with all bond interactions with the class method -single() implemented and all angle interactions with the class method -born_matrix() implemented. +improper intramolecular interactions. Related commands """""""""""""""" From d84ee0c4f1f97728d69a55163d8356eba16e290b Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 28 Sep 2023 17:48:17 +0300 Subject: [PATCH 42/60] Update compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 223e4da66a..41b5f64a67 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -325,6 +325,9 @@ void ComputeStressMopProfile::compute_array() } } + // sum dihedral contribution over all procs + MPI_Allreduce(&dihedral_local[0][0],&dihedral_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + for (int ibin = 0; ibin < nbins; ibin++) { array[ibin][0] = coord[ibin]; From 5830dec742dc5c847a72b4bb1e6c7b5fb2aea756 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 16 Oct 2023 15:09:10 -0600 Subject: [PATCH 43/60] new compute: reaxff/bonds/local --- src/REAXFF/compute_reaxff_bonds_local.cpp | 217 ++++++++++++++++++++++ src/REAXFF/compute_reaxff_bonds_local.h | 60 ++++++ 2 files changed, 277 insertions(+) create mode 100644 src/REAXFF/compute_reaxff_bonds_local.cpp create mode 100644 src/REAXFF/compute_reaxff_bonds_local.h diff --git a/src/REAXFF/compute_reaxff_bonds_local.cpp b/src/REAXFF/compute_reaxff_bonds_local.cpp new file mode 100644 index 0000000000..c32b404e86 --- /dev/null +++ b/src/REAXFF/compute_reaxff_bonds_local.cpp @@ -0,0 +1,217 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#include "compute_reaxff_bonds_local.h" +#include "atom.h" +#include "molecule.h" +#include "update.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "neigh_list.h" + +#include "pair_reaxff.h" +#include "reaxff_api.h" + +using namespace LAMMPS_NS; +using namespace ReaxFF; + +/* ---------------------------------------------------------------------- */ + +ComputeReaxFFBondsLocal::ComputeReaxFFBondsLocal(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + alocal(nullptr), abo(nullptr), neighid(nullptr), numneigh(nullptr), reaxff(nullptr) +{ + if (atom->tag_consecutive() == 0) + error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds/local"); + + local_flag = 1; + + nvalues = 7 + 2*MAXREAXBOND; + prev_nvalues = 0; + + // initialize output + + nlocal = atom->nlocal; + + size_local_rows = atom->nlocal; + size_local_cols = 7 + 2*MAXREAXBOND; + + allocate(nlocal); +} + + +/* ---------------------------------------------------------------------- */ + +ComputeReaxFFBondsLocal::~ComputeReaxFFBondsLocal() +{ + memory->destroy(alocal); + destroy(); +} + +void ComputeReaxFFBondsLocal::destroy() +{ + memory->destroy(abo); + memory->destroy(neighid); + memory->destroy(numneigh); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFBondsLocal::allocate(int n) +{ + memory->create(abo,n,MAXREAXBOND,"reaxff/bonds/local:abo"); + memory->create(neighid,n,MAXREAXBOND,"reaxff/bonds/local:neighid"); + memory->create(numneigh,n,"reaxff/bonds/local:numneigh"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFBondsLocal::init() +{ + reaxff = dynamic_cast(force->pair_match("^reax..",0)); + if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/bonds/local without " + "pair_style reaxff, reaxff/kk, or reaxff/omp"); +} + +/* ---------------------------------------------------------------------- */ + +int ComputeReaxFFBondsLocal::FindBond() +{ + int *ilist, i, ii, inum; + int j, pj, nj; + tagint jtag; + double bo_tmp,bo_cut; + + inum = reaxff->list->inum; + ilist = reaxff->list->ilist; + bond_data *bo_ij; + bo_cut = reaxff->api->control->bg_cut; + + tagint *tag = atom->tag; + int numbonds = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + nj = 0; + + for (pj = Start_Index(i, reaxff->api->lists); pj < End_Index(i, reaxff->api->lists); ++pj) { + bo_ij = &(reaxff->api->lists->select.bond_list[pj]); + j = bo_ij->nbr; + jtag = tag[j]; + bo_tmp = bo_ij->bo_data.BO; + + if (bo_tmp > bo_cut) { + neighid[i][nj] = jtag; + abo[i][nj] = bo_tmp; + nj ++; + } + } + numneigh[i] = nj; + if (nj > numbonds) numbonds = nj; + } + return numbonds; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFBondsLocal::compute_local() +{ + invoked_local = update->ntimestep; + + // count local entries and compute bond info + if (atom->nlocal > nlocal) { + destroy(); + allocate(atom->nlocal); + } + + { + const int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + numneigh[i] = 0; + for (int j = 0; j < MAXREAXBOND; j++) { + neighid[i][j] = 0; + abo[i][j] = 0.0; + } + } + } + + int maxnumbonds = FindBond(); + nvalues = 7+2*maxnumbonds; + + if(atom->nlocal > nlocal || nvalues > prev_nvalues) { + reallocate(); + } + + size_local_rows = nlocal; + size_local_cols = nvalues; + + for (int i = 0; i < nlocal; ++i) { + auto ptr = alocal[i]; + int numbonds = numneigh[i]; + ptr[0] = atom->tag[i]; + ptr[1] = atom->type[i]; + ptr[2] = numbonds; + + int j = 3; + + for (int k = 0; k < numbonds; k++) { + ptr[j++] = neighid[i][k]; + } + + ptr[j++] = (atom->molecule == nullptr) ? 0.0 : atom->molecule[i]; + + for (int k = 0; k < numbonds; k++) { + ptr[j++] = abo[i][k]; + } + + ptr[j++] = reaxff->api->workspace->total_bond_order[i]; + ptr[j++] = reaxff->api->workspace->nlp[i]; + ptr[j++] = atom->q[i]; + + // clear any remaining + for(; j < nvalues; ++j) { + ptr[j] = 0.0; + } + } +} + +void ComputeReaxFFBondsLocal::reallocate() +{ + nlocal = atom->nlocal; + + // grow array_local + memory->destroy(alocal); + memory->create(alocal,nlocal,nvalues,"reaxff/bonds/local:array_local"); + array_local = alocal; + + prev_nvalues = nvalues; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeReaxFFBondsLocal::memory_usage() +{ + double bytes = (double)nlocal*nvalues * sizeof(double); + bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); + bytes += (double)(nlocal) * sizeof(int); + return bytes; +} diff --git a/src/REAXFF/compute_reaxff_bonds_local.h b/src/REAXFF/compute_reaxff_bonds_local.h new file mode 100644 index 0000000000..2347db3a26 --- /dev/null +++ b/src/REAXFF/compute_reaxff_bonds_local.h @@ -0,0 +1,60 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(reaxff/bonds/local,ComputeReaxFFBondsLocal); +// clang-format on +#else + +#ifndef LMP_COMPUTE_REAXFF_BONDS_LOCAL_H +#define LMP_COMPUTE_REAXFF_BONDS_LOCAL_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeReaxFFBondsLocal : public Compute { + public: + ComputeReaxFFBondsLocal(class LAMMPS *, int, char **); + ~ComputeReaxFFBondsLocal() override; + void init() override; + void compute_local() override; + double memory_usage() override; + + private: + int nlocal; + int nvalues; + int prev_nvalues; + + double **alocal; + tagint **neighid; + double **abo; + int *numneigh; + class PairReaxFF *reaxff; + + int FindBond(); + + void allocate(int); + void destroy(); + void reallocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif From bbc2794df29920bc1673ac0501518669fcfb822d Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 16 Oct 2023 15:09:16 -0600 Subject: [PATCH 44/60] add reaxff/bonds/local/kk --- .../compute_reaxff_bonds_local_kokkos.cpp | 125 ++++++++++++++++++ .../compute_reaxff_bonds_local_kokkos.h | 68 ++++++++++ src/KOKKOS/pair_reaxff_kokkos.cpp | 55 ++++++++ src/KOKKOS/pair_reaxff_kokkos.h | 19 +++ 4 files changed, 267 insertions(+) create mode 100644 src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp create mode 100644 src/KOKKOS/compute_reaxff_bonds_local_kokkos.h diff --git a/src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp b/src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp new file mode 100644 index 0000000000..a56d243ed6 --- /dev/null +++ b/src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp @@ -0,0 +1,125 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#include "compute_reaxff_bonds_local_kokkos.h" +#include "atom.h" +#include "molecule.h" +#include "update.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "neigh_list.h" + +#include "memory_kokkos.h" +#include "pair_reaxff_kokkos.h" +#include "reaxff_api.h" + +using namespace LAMMPS_NS; +using namespace ReaxFF; + +/* ---------------------------------------------------------------------- */ + +template +ComputeReaxFFBondsLocalKokkos::ComputeReaxFFBondsLocalKokkos(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + alocal(nullptr), reaxff(nullptr) +{ + if (atom->tag_consecutive() == 0) + error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds/local"); + + local_flag = 1; + + nvalues = 7 + 2*MAXREAXBOND; + prev_nvalues = 0; + + // initialize output + + nlocal = atom->nlocal; + + size_local_rows = atom->nlocal; + size_local_cols = 7 + 2*MAXREAXBOND; + printf("RUNNING KOKKOS VERSION\n"); +} + + +/* ---------------------------------------------------------------------- */ + +template +ComputeReaxFFBondsLocalKokkos::~ComputeReaxFFBondsLocalKokkos() +{ + memoryKK->destroy_kokkos(k_alocal, alocal); +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFBondsLocalKokkos::init() +{ + reaxff = dynamic_cast(force->pair_match("^reax../kk",0)); + if (reaxff == nullptr) error->all(FLERR,"Cannot use fix reaxff/bonds without " + "pair_style reaxff/kk"); +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFBondsLocalKokkos::compute_local() +{ + invoked_local = update->ntimestep; + + int maxnumbonds = 0; + if (reaxff->execution_space == Device) + device_pair()->FindBond(maxnumbonds); + else + host_pair()->FindBond(maxnumbonds); + nvalues = 7+2*maxnumbonds; + + if(atom->nlocal > nlocal || nvalues > prev_nvalues) { + nlocal = atom->nlocal; + memoryKK->destroy_kokkos(k_alocal, alocal); + memoryKK->create_kokkos(k_alocal, alocal, atom->nlocal, nvalues,"reaxff/bonds/local:alocal"); + prev_nvalues = nvalues; + array_local = alocal; + } + + size_local_rows = nlocal; + size_local_cols = nvalues; + + if (reaxff->execution_space == Device) + device_pair()->PackBondInfo(k_alocal); + else + host_pair()->PackBondInfo(k_alocal); +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +template +double ComputeReaxFFBondsLocalKokkos::memory_usage() +{ + double bytes = (double)nlocal*nvalues * sizeof(double); + return bytes; +} + +namespace LAMMPS_NS { +template class ComputeReaxFFBondsLocalKokkos; +#ifdef LMP_KOKKOS_GPU +template class ComputeReaxFFBondsLocalKokkos; +#endif +} diff --git a/src/KOKKOS/compute_reaxff_bonds_local_kokkos.h b/src/KOKKOS/compute_reaxff_bonds_local_kokkos.h new file mode 100644 index 0000000000..bfb5521199 --- /dev/null +++ b/src/KOKKOS/compute_reaxff_bonds_local_kokkos.h @@ -0,0 +1,68 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(reaxff/bonds/local/kk,ComputeReaxFFBondsLocalKokkos); +ComputeStyle(reaxff/bonds/local/kk/device,ComputeReaxFFBondsLocalKokkos); +ComputeStyle(reaxff/bonds/local/kk/host,ComputeReaxFFBondsLocalKokkos); +// clang-format on +#else + +#ifndef LMP_COMPUTE_REAXFF_BONDS_LOCAL_KOKKOS_H +#define LMP_COMPUTE_REAXFF_BONDS_LOCAL_KOKKOS_H + +#include "compute_reaxff_bonds_local.h" +#include "pair_reaxff_kokkos.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +template +class ComputeReaxFFBondsLocalKokkos : public Compute { + public: + using device_type = DeviceType; + using AT = ArrayTypes; + + ComputeReaxFFBondsLocalKokkos(class LAMMPS *, int, char **); + ~ComputeReaxFFBondsLocalKokkos() override; + void init() override; + void compute_local() override; + double memory_usage() override; + + private: + int nlocal; + int nvalues; + int prev_nvalues; + + double **alocal; + typename AT::tdual_float_2d k_alocal; + PairReaxFF *reaxff; + + auto device_pair() { + return dynamic_cast*>(reaxff); + } + + auto host_pair() { + return dynamic_cast*>(reaxff); + } +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index c7d54b80cd..fbf90f140d 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4290,6 +4290,61 @@ void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bo /* ---------------------------------------------------------------------- */ +template +void PairReaxFFKokkos::PackBondInfo(DAT::tdual_float_2d k_alocal) +{ + d_alocal = k_alocal.view(); + k_params_sing.template sync(); + atomKK->sync(execution_space,TAG_MASK|TYPE_MASK|Q_MASK|MOLECULE_MASK); + + tag = atomKK->k_tag.view(); + type = atomKK->k_type.view(); + q = atomKK->k_q.view(); + if (atom->molecule) + molecule = atomKK->k_molecule.view(); + + copymode = 1; + nlocal = atomKK->nlocal; + PairReaxKokkosPackBondInfoFunctor pack_bond_info_functor(this); + Kokkos::parallel_for(nlocal, pack_bond_info_functor); + copymode = 0; + + k_alocal.modify(); + k_alocal.sync(); +} + +template +KOKKOS_INLINE_FUNCTION +void PairReaxFFKokkos::pack_bond_info_item(const int i) const +{ + const int numbonds = d_numneigh_bonds[i]; + d_alocal(i,0) = tag[i]; + d_alocal(i,1) = type[i]; + d_alocal(i,2) = numbonds; + + int j = 3; + + for (int k = 0; k < numbonds; k++) { + d_alocal(i,j++) = d_neighid(i,k); + } + + d_alocal(i,j++) = molecule.data() ? molecule[i] : 0.0; + + for (int k = 0; k < numbonds; k++) { + d_alocal(i,j++) = d_abo(i,k); + } + + d_alocal(i,j++) = d_total_bo[i]; + d_alocal(i,j++) = paramssing(type[i]).nlp_opt - d_Delta_lp[i]; + d_alocal(i,j++) = q[i]; + + for(; j < d_alocal.extent(1); ++j) { + d_alocal(i,j) = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + template void PairReaxFFKokkos::FindBondSpecies() { diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index 1ad0955a1e..32007e9970 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -135,6 +135,7 @@ class PairReaxFFKokkos : public PairReaxFF { double memory_usage(); void FindBond(int &); void PackBondBuffer(DAT::tdual_ffloat_1d, int &); + void PackBondInfo(DAT::tdual_float_2d); void FindBondSpecies(); template @@ -292,6 +293,9 @@ class PairReaxFFKokkos : public PairReaxFF { KOKKOS_INLINE_FUNCTION void pack_bond_buffer_item(int, int&, const bool&) const; + KOKKOS_INLINE_FUNCTION + void pack_bond_info_item(const int) const; + KOKKOS_INLINE_FUNCTION void operator()(TagPairReaxFindBondSpeciesZero, const int&) const; @@ -506,6 +510,8 @@ class PairReaxFFKokkos : public PairReaxFF { typename AT::t_ffloat_1d d_buf; DAT::tdual_int_scalar k_nbuf_local; + typename AT::t_float_2d d_alocal; + typedef Kokkos::View t_reax_int4_2d; t_reax_int4_2d d_angular_pack, d_torsion_pack; @@ -549,6 +555,19 @@ struct PairReaxKokkosPackBondBufferFunctor { } }; +template +struct PairReaxKokkosPackBondInfoFunctor { + using device_type = DeviceType; + using value_type = int; + PairReaxFFKokkos c; + PairReaxKokkosPackBondInfoFunctor(PairReaxFFKokkos* c_ptr):c(*c_ptr) {}; + + KOKKOS_INLINE_FUNCTION + void operator()(const int i) const { + c.pack_bond_info_item(i); + } +}; + } #endif From fea5f5a2438d5ef1c563e945783f93a9560a18f4 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 16 Oct 2023 16:52:34 -0600 Subject: [PATCH 45/60] First serial version of Steve's suggestion --- examples/reaxff/in.reaxff.tatb | 9 +- ...nds_local.cpp => compute_reaxff_bonds.cpp} | 174 +++++++++--------- ...f_bonds_local.h => compute_reaxff_bonds.h} | 27 ++- 3 files changed, 107 insertions(+), 103 deletions(-) rename src/REAXFF/{compute_reaxff_bonds_local.cpp => compute_reaxff_bonds.cpp} (57%) rename src/REAXFF/{compute_reaxff_bonds_local.h => compute_reaxff_bonds.h} (73%) diff --git a/examples/reaxff/in.reaxff.tatb b/examples/reaxff/in.reaxff.tatb index 6cf7828cf1..fd0d846154 100644 --- a/examples/reaxff/in.reaxff.tatb +++ b/examples/reaxff/in.reaxff.tatb @@ -31,8 +31,15 @@ neigh_modify delay 0 every 5 check no fix 1 all nve fix 2 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff fix 4 all reaxff/bonds 5 bonds.reaxff +compute bonds all reaxff/bonds variable nqeq equal f_2 +# dumps out the local bond information +dump 1 all local 5 bonds_compute.reaxff c_bonds[1] c_bonds[2] c_bonds[3] + +# dumps out the peratom bond information +dump 2 all custom 5 bonds_atom.reaxff c_bonds[*] + thermo 5 thermo_style custom step temp epair etotal press & v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa & @@ -50,6 +57,6 @@ timestep 0.0625 # axes yes 0.8 0.02 view 60 -30 #dump_modify 3 pad 3 -fix 3 all reaxff/species 1 5 5 species.tatb +#fix 3 all reaxff/species 1 5 5 species.tatb run 25 diff --git a/src/REAXFF/compute_reaxff_bonds_local.cpp b/src/REAXFF/compute_reaxff_bonds.cpp similarity index 57% rename from src/REAXFF/compute_reaxff_bonds_local.cpp rename to src/REAXFF/compute_reaxff_bonds.cpp index c32b404e86..7f1605263c 100644 --- a/src/REAXFF/compute_reaxff_bonds_local.cpp +++ b/src/REAXFF/compute_reaxff_bonds.cpp @@ -16,7 +16,7 @@ Contributing author: Richard Berger (LANL) ------------------------------------------------------------------------- */ -#include "compute_reaxff_bonds_local.h" +#include "compute_reaxff_bonds.h" #include "atom.h" #include "molecule.h" #include "update.h" @@ -33,39 +33,36 @@ using namespace ReaxFF; /* ---------------------------------------------------------------------- */ -ComputeReaxFFBondsLocal::ComputeReaxFFBondsLocal(LAMMPS *lmp, int narg, char **arg) : +ComputeReaxFFBonds::ComputeReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - alocal(nullptr), abo(nullptr), neighid(nullptr), numneigh(nullptr), reaxff(nullptr) + abo(nullptr), neighid(nullptr), numneigh(nullptr), reaxff(nullptr) { if (atom->tag_consecutive() == 0) - error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds/local"); + error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds"); local_flag = 1; - - nvalues = 7 + 2*MAXREAXBOND; - prev_nvalues = 0; + peratom_flag = 1; // initialize output - nlocal = atom->nlocal; + nlocal = -1; + prev_nbonds = -1; - size_local_rows = atom->nlocal; - size_local_cols = 7 + 2*MAXREAXBOND; + size_peratom_cols = 7; - allocate(nlocal); + size_local_rows = 0; + size_local_cols = 3; + + invoked_bonds = -1; } /* ---------------------------------------------------------------------- */ -ComputeReaxFFBondsLocal::~ComputeReaxFFBondsLocal() -{ - memory->destroy(alocal); - destroy(); -} - -void ComputeReaxFFBondsLocal::destroy() +ComputeReaxFFBonds::~ComputeReaxFFBonds() { + memory->destroy(array_local); + memory->destroy(array_atom); memory->destroy(abo); memory->destroy(neighid); memory->destroy(numneigh); @@ -73,25 +70,16 @@ void ComputeReaxFFBondsLocal::destroy() /* ---------------------------------------------------------------------- */ -void ComputeReaxFFBondsLocal::allocate(int n) -{ - memory->create(abo,n,MAXREAXBOND,"reaxff/bonds/local:abo"); - memory->create(neighid,n,MAXREAXBOND,"reaxff/bonds/local:neighid"); - memory->create(numneigh,n,"reaxff/bonds/local:numneigh"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeReaxFFBondsLocal::init() +void ComputeReaxFFBonds::init() { reaxff = dynamic_cast(force->pair_match("^reax..",0)); - if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/bonds/local without " - "pair_style reaxff, reaxff/kk, or reaxff/omp"); + if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/bonds without " + "pair_style reaxff, reaxff/kk, or reaxff/omp"); } /* ---------------------------------------------------------------------- */ -int ComputeReaxFFBondsLocal::FindBond() +int ComputeReaxFFBonds::FindBond() { int *ilist, i, ii, inum; int j, pj, nj; @@ -119,98 +107,108 @@ int ComputeReaxFFBondsLocal::FindBond() if (bo_tmp > bo_cut) { neighid[i][nj] = jtag; abo[i][nj] = bo_tmp; - nj ++; + nj++; } } numneigh[i] = nj; - if (nj > numbonds) numbonds = nj; + numbonds += nj; } return numbonds; } + /* ---------------------------------------------------------------------- */ -void ComputeReaxFFBondsLocal::compute_local() +void ComputeReaxFFBonds::compute_bonds() +{ + invoked_bonds = update->ntimestep; + + if (atom->nlocal > nlocal) { + memory->destroy(abo); + memory->destroy(neighid); + memory->destroy(numneigh); + memory->destroy(array_atom); + nlocal = atom->nlocal; + memory->create(abo, nlocal, MAXREAXBOND, "reaxff/bonds:abo"); + memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/bonds:neighid"); + memory->create(numneigh, nlocal, "reaxff/bonds:numneigh"); + memory->create(array_atom, nlocal, 7, "reaxff/bonds:array_atom"); + } + + for (int i = 0; i < nlocal; i++) { + numneigh[i] = 0; + for (int j = 0; j < MAXREAXBOND; j++) { + neighid[i][j] = 0; + abo[i][j] = 0.0; + } + } + + nbonds = FindBond(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFBonds::compute_local() { invoked_local = update->ntimestep; - // count local entries and compute bond info - if (atom->nlocal > nlocal) { - destroy(); - allocate(atom->nlocal); + if(invoked_bonds < update->ntimestep) { + compute_bonds(); } - { - const int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) { - numneigh[i] = 0; - for (int j = 0; j < MAXREAXBOND; j++) { - neighid[i][j] = 0; - abo[i][j] = 0.0; - } - } + if(nbonds > prev_nbonds) { + // grow array_local + memory->destroy(array_local); + memory->create(array_local, nbonds, 3, "reaxff/bonds:array_local"); + prev_nbonds = nbonds; } - int maxnumbonds = FindBond(); - nvalues = 7+2*maxnumbonds; + size_local_rows = nbonds; - if(atom->nlocal > nlocal || nvalues > prev_nvalues) { - reallocate(); - } - - size_local_rows = nlocal; - size_local_cols = nvalues; + int b = 0; for (int i = 0; i < nlocal; ++i) { - auto ptr = alocal[i]; - int numbonds = numneigh[i]; - ptr[0] = atom->tag[i]; - ptr[1] = atom->type[i]; - ptr[2] = numbonds; - - int j = 3; + const int numbonds = numneigh[i]; for (int k = 0; k < numbonds; k++) { - ptr[j++] = neighid[i][k]; - } - - ptr[j++] = (atom->molecule == nullptr) ? 0.0 : atom->molecule[i]; - - for (int k = 0; k < numbonds; k++) { - ptr[j++] = abo[i][k]; - } - - ptr[j++] = reaxff->api->workspace->total_bond_order[i]; - ptr[j++] = reaxff->api->workspace->nlp[i]; - ptr[j++] = atom->q[i]; - - // clear any remaining - for(; j < nvalues; ++j) { - ptr[j] = 0.0; + auto bond = array_local[b++]; + bond[0] = i; + bond[1] = neighid[i][k]; + bond[2] = abo[i][k]; } } } -void ComputeReaxFFBondsLocal::reallocate() +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFBonds::compute_peratom() { - nlocal = atom->nlocal; + invoked_peratom = update->ntimestep; - // grow array_local - memory->destroy(alocal); - memory->create(alocal,nlocal,nvalues,"reaxff/bonds/local:array_local"); - array_local = alocal; + if(invoked_bonds < update->ntimestep) { + compute_bonds(); + } - prev_nvalues = nvalues; + for (int i = 0; i < nlocal; ++i) { + auto ptr = array_atom[i]; + ptr[0] = atom->tag[i]; + ptr[1] = atom->type[i]; + ptr[2] = numneigh[i]; + ptr[3] = (atom->molecule == nullptr) ? 0.0 : atom->molecule[i]; + ptr[4] = reaxff->api->workspace->total_bond_order[i]; + ptr[5] = reaxff->api->workspace->nlp[i]; + ptr[6] = atom->q[i]; + } } /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ -double ComputeReaxFFBondsLocal::memory_usage() +double ComputeReaxFFBonds::memory_usage() { - double bytes = (double)nlocal*nvalues * sizeof(double); + double bytes = (double)(nbonds*3) * sizeof(double); + bytes += (double)(nlocal*7) * sizeof(double); bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); bytes += (double)(nlocal) * sizeof(int); return bytes; diff --git a/src/REAXFF/compute_reaxff_bonds_local.h b/src/REAXFF/compute_reaxff_bonds.h similarity index 73% rename from src/REAXFF/compute_reaxff_bonds_local.h rename to src/REAXFF/compute_reaxff_bonds.h index 2347db3a26..6b00cef7ed 100644 --- a/src/REAXFF/compute_reaxff_bonds_local.h +++ b/src/REAXFF/compute_reaxff_bonds.h @@ -17,41 +17,40 @@ #ifdef COMPUTE_CLASS // clang-format off -ComputeStyle(reaxff/bonds/local,ComputeReaxFFBondsLocal); +ComputeStyle(reaxff/bonds,ComputeReaxFFBonds); // clang-format on #else -#ifndef LMP_COMPUTE_REAXFF_BONDS_LOCAL_H -#define LMP_COMPUTE_REAXFF_BONDS_LOCAL_H +#ifndef LMP_COMPUTE_REAXFF_BONDS_H +#define LMP_COMPUTE_REAXFF_BONDS_H #include "compute.h" namespace LAMMPS_NS { -class ComputeReaxFFBondsLocal : public Compute { +class ComputeReaxFFBonds : public Compute { public: - ComputeReaxFFBondsLocal(class LAMMPS *, int, char **); - ~ComputeReaxFFBondsLocal() override; + ComputeReaxFFBonds(class LAMMPS *, int, char **); + ~ComputeReaxFFBonds() override; void init() override; void compute_local() override; + void compute_peratom() override; + virtual void compute_bonds(); double memory_usage() override; - private: + protected: + bigint invoked_bonds; // last timestep on which compute_bonds() was invoked int nlocal; - int nvalues; - int prev_nvalues; + int nbonds; + int prev_nbonds; - double **alocal; tagint **neighid; double **abo; int *numneigh; class PairReaxFF *reaxff; + private: int FindBond(); - - void allocate(int); - void destroy(); - void reallocate(); }; } // namespace LAMMPS_NS From a72a3ed50d12365ebb41a02001454fe667036b23 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Tue, 17 Oct 2023 16:24:14 -0600 Subject: [PATCH 46/60] Kokkos version of compute reaxff/bonds --- src/KOKKOS/compute_reaxff_bonds_kokkos.cpp | 192 ++++++++++++++++++ ...kokkos.h => compute_reaxff_bonds_kokkos.h} | 30 ++- .../compute_reaxff_bonds_local_kokkos.cpp | 125 ------------ src/KOKKOS/pair_reaxff_kokkos.cpp | 55 ----- src/KOKKOS/pair_reaxff_kokkos.h | 19 -- 5 files changed, 206 insertions(+), 215 deletions(-) create mode 100644 src/KOKKOS/compute_reaxff_bonds_kokkos.cpp rename src/KOKKOS/{compute_reaxff_bonds_local_kokkos.h => compute_reaxff_bonds_kokkos.h} (67%) delete mode 100644 src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp diff --git a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp b/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp new file mode 100644 index 0000000000..921acb9193 --- /dev/null +++ b/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp @@ -0,0 +1,192 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#include "compute_reaxff_bonds_kokkos.h" +#include "atom.h" +#include "molecule.h" +#include "update.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "neigh_list.h" + +#include "memory_kokkos.h" +#include "pair_reaxff_kokkos.h" +#include "reaxff_api.h" + +using namespace LAMMPS_NS; +using namespace ReaxFF; + +/* ---------------------------------------------------------------------- */ + +template +ComputeReaxFFBondsKokkos::ComputeReaxFFBondsKokkos(LAMMPS *lmp, int narg, char **arg) : + ComputeReaxFFBonds(lmp, narg, arg), + nbuf(-1), buf(nullptr) +{ +} + + +/* ---------------------------------------------------------------------- */ + +template +ComputeReaxFFBondsKokkos::~ComputeReaxFFBondsKokkos() +{ + memoryKK->destroy_kokkos(k_buf, buf); +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFBondsKokkos::init() +{ + reaxff = dynamic_cast(force->pair_match("^reax../kk",0)); + if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/bonds without " + "pair_style reaxff/kk"); +} + +/* ---------------------------------------------------------------------- */ +template +void ComputeReaxFFBondsKokkos::compute_bonds() +{ + if (atom->nlocal > nlocal) { + memory->destroy(array_atom); + nlocal = atom->nlocal; + memory->create(array_atom, nlocal, 7, "reaxff/bonds:array_atom"); + } + + // retrieve bond information from kokkos pair style. the data potentially + // lives on device. it is copied into buf on the host in a condensed format + // compute_local and compute_atom then expand the data from this buffer into + // appropiate arrays for consumption by others (e.g. dump local, dump custom + // or library interface) + int maxnumbonds = 0; + if (reaxff->execution_space == Device) + device_pair()->FindBond(maxnumbonds); + else + host_pair()->FindBond(maxnumbonds); + + nbuf = 1+(maxnumbonds*2 + 7)*nlocal; + + if(!buf || k_buf.extent(0) < nbuf) { + memoryKK->destroy_kokkos(k_buf, buf); + memoryKK->create_kokkos(k_buf, buf, nbuf, "reaxff/bonds:buf"); + } + + // Pass information to buffer, will sync to host + int nbuf_local; + if (reaxff->execution_space == Device) + device_pair()->PackBondBuffer(k_buf, nbuf_local); + else + host_pair()->PackBondBuffer(k_buf, nbuf_local); + buf[0] = nlocal; + + // Extract number of bonds from buffer + nbonds = 0; + int j = 1; + for (int i = 0; i < nlocal; i++) { + int numbonds = static_cast(buf[j+5]); + nbonds += numbonds; + j += 2*numbonds + 7; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFBondsKokkos::compute_local() +{ + invoked_local = update->ntimestep; + + if(invoked_bonds < update->ntimestep) { + compute_bonds(); + } + + if(nbonds > prev_nbonds) { + // grow array_local + memory->destroy(array_local); + memory->create(array_local, nbonds, 3, "reaxff/bonds:array_local"); + prev_nbonds = nbonds; + } + + size_local_rows = nbonds; + + // extract local bond information from buffer + int b = 0; + int j = 1; + + for (int i = 0; i < nlocal; ++i) { + const int numbonds = static_cast(buf[j+5]); + const int neigh_offset = j + 6; + const int bo_offset = neigh_offset + numbonds + 1; + for (int k = 0; k < numbonds; k++) { + auto bond = array_local[b++]; + bond[0] = i; + bond[1] = static_cast (buf[neigh_offset+k]); + bond[2] = buf[bo_offset+k]; + } + j += 2*numbonds + 7; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFBondsKokkos::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + if(invoked_bonds < update->ntimestep) { + compute_bonds(); + } + + // extract peratom bond information from buffer + int j = 1; + for (int i = 0; i < nlocal; ++i) { + auto ptr = array_atom[i]; + int numbonds = static_cast(buf[j+5]); + const int mol_offset = j + 6 + numbonds; + ptr[0] = buf[j]; // jtag + ptr[1] = buf[j+1]; // itype + ptr[2] = numbonds; + ptr[3] = buf[mol_offset]; // mol + ptr[4] = buf[j+2]; // sbo + ptr[5] = buf[j+3]; // nlp + ptr[6] = buf[j+4]; // q + j += 2*numbonds + 7; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +template +double ComputeReaxFFBondsKokkos::memory_usage() +{ + double bytes = (double)(nbonds*3) * sizeof(double); + bytes += (double)(nlocal*7) * sizeof(double); + return bytes; +} + +namespace LAMMPS_NS { +template class ComputeReaxFFBondsKokkos; +#ifdef LMP_KOKKOS_GPU +template class ComputeReaxFFBondsKokkos; +#endif +} diff --git a/src/KOKKOS/compute_reaxff_bonds_local_kokkos.h b/src/KOKKOS/compute_reaxff_bonds_kokkos.h similarity index 67% rename from src/KOKKOS/compute_reaxff_bonds_local_kokkos.h rename to src/KOKKOS/compute_reaxff_bonds_kokkos.h index bfb5521199..45020ffa81 100644 --- a/src/KOKKOS/compute_reaxff_bonds_local_kokkos.h +++ b/src/KOKKOS/compute_reaxff_bonds_kokkos.h @@ -17,41 +17,39 @@ #ifdef COMPUTE_CLASS // clang-format off -ComputeStyle(reaxff/bonds/local/kk,ComputeReaxFFBondsLocalKokkos); -ComputeStyle(reaxff/bonds/local/kk/device,ComputeReaxFFBondsLocalKokkos); -ComputeStyle(reaxff/bonds/local/kk/host,ComputeReaxFFBondsLocalKokkos); +ComputeStyle(reaxff/bonds/kk,ComputeReaxFFBondsKokkos); +ComputeStyle(reaxff/bonds/kk/device,ComputeReaxFFBondsKokkos); +ComputeStyle(reaxff/bonds/kk/host,ComputeReaxFFBondsKokkos); // clang-format on #else -#ifndef LMP_COMPUTE_REAXFF_BONDS_LOCAL_KOKKOS_H -#define LMP_COMPUTE_REAXFF_BONDS_LOCAL_KOKKOS_H +#ifndef LMP_COMPUTE_REAXFF_BONDS_KOKKOS_H +#define LMP_COMPUTE_REAXFF_BONDS_KOKKOS_H -#include "compute_reaxff_bonds_local.h" +#include "compute_reaxff_bonds.h" #include "pair_reaxff_kokkos.h" #include "kokkos_type.h" namespace LAMMPS_NS { template -class ComputeReaxFFBondsLocalKokkos : public Compute { +class ComputeReaxFFBondsKokkos : public ComputeReaxFFBonds { public: using device_type = DeviceType; using AT = ArrayTypes; - ComputeReaxFFBondsLocalKokkos(class LAMMPS *, int, char **); - ~ComputeReaxFFBondsLocalKokkos() override; + ComputeReaxFFBondsKokkos(class LAMMPS *, int, char **); + ~ComputeReaxFFBondsKokkos() override; void init() override; void compute_local() override; + void compute_peratom() override; + void compute_bonds() override; double memory_usage() override; private: - int nlocal; - int nvalues; - int prev_nvalues; - - double **alocal; - typename AT::tdual_float_2d k_alocal; - PairReaxFF *reaxff; + int nbuf; + double *buf; + typename AT::tdual_float_1d k_buf; auto device_pair() { return dynamic_cast*>(reaxff); diff --git a/src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp b/src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp deleted file mode 100644 index a56d243ed6..0000000000 --- a/src/KOKKOS/compute_reaxff_bonds_local_kokkos.cpp +++ /dev/null @@ -1,125 +0,0 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - LAMMPS development team: developers@lammps.org - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Richard Berger (LANL) -------------------------------------------------------------------------- */ - -#include "compute_reaxff_bonds_local_kokkos.h" -#include "atom.h" -#include "molecule.h" -#include "update.h" -#include "force.h" -#include "memory.h" -#include "error.h" -#include "neigh_list.h" - -#include "memory_kokkos.h" -#include "pair_reaxff_kokkos.h" -#include "reaxff_api.h" - -using namespace LAMMPS_NS; -using namespace ReaxFF; - -/* ---------------------------------------------------------------------- */ - -template -ComputeReaxFFBondsLocalKokkos::ComputeReaxFFBondsLocalKokkos(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - alocal(nullptr), reaxff(nullptr) -{ - if (atom->tag_consecutive() == 0) - error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds/local"); - - local_flag = 1; - - nvalues = 7 + 2*MAXREAXBOND; - prev_nvalues = 0; - - // initialize output - - nlocal = atom->nlocal; - - size_local_rows = atom->nlocal; - size_local_cols = 7 + 2*MAXREAXBOND; - printf("RUNNING KOKKOS VERSION\n"); -} - - -/* ---------------------------------------------------------------------- */ - -template -ComputeReaxFFBondsLocalKokkos::~ComputeReaxFFBondsLocalKokkos() -{ - memoryKK->destroy_kokkos(k_alocal, alocal); -} - -/* ---------------------------------------------------------------------- */ - -template -void ComputeReaxFFBondsLocalKokkos::init() -{ - reaxff = dynamic_cast(force->pair_match("^reax../kk",0)); - if (reaxff == nullptr) error->all(FLERR,"Cannot use fix reaxff/bonds without " - "pair_style reaxff/kk"); -} - -/* ---------------------------------------------------------------------- */ - -template -void ComputeReaxFFBondsLocalKokkos::compute_local() -{ - invoked_local = update->ntimestep; - - int maxnumbonds = 0; - if (reaxff->execution_space == Device) - device_pair()->FindBond(maxnumbonds); - else - host_pair()->FindBond(maxnumbonds); - nvalues = 7+2*maxnumbonds; - - if(atom->nlocal > nlocal || nvalues > prev_nvalues) { - nlocal = atom->nlocal; - memoryKK->destroy_kokkos(k_alocal, alocal); - memoryKK->create_kokkos(k_alocal, alocal, atom->nlocal, nvalues,"reaxff/bonds/local:alocal"); - prev_nvalues = nvalues; - array_local = alocal; - } - - size_local_rows = nlocal; - size_local_cols = nvalues; - - if (reaxff->execution_space == Device) - device_pair()->PackBondInfo(k_alocal); - else - host_pair()->PackBondInfo(k_alocal); -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -template -double ComputeReaxFFBondsLocalKokkos::memory_usage() -{ - double bytes = (double)nlocal*nvalues * sizeof(double); - return bytes; -} - -namespace LAMMPS_NS { -template class ComputeReaxFFBondsLocalKokkos; -#ifdef LMP_KOKKOS_GPU -template class ComputeReaxFFBondsLocalKokkos; -#endif -} diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index fbf90f140d..c7d54b80cd 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4290,61 +4290,6 @@ void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bo /* ---------------------------------------------------------------------- */ -template -void PairReaxFFKokkos::PackBondInfo(DAT::tdual_float_2d k_alocal) -{ - d_alocal = k_alocal.view(); - k_params_sing.template sync(); - atomKK->sync(execution_space,TAG_MASK|TYPE_MASK|Q_MASK|MOLECULE_MASK); - - tag = atomKK->k_tag.view(); - type = atomKK->k_type.view(); - q = atomKK->k_q.view(); - if (atom->molecule) - molecule = atomKK->k_molecule.view(); - - copymode = 1; - nlocal = atomKK->nlocal; - PairReaxKokkosPackBondInfoFunctor pack_bond_info_functor(this); - Kokkos::parallel_for(nlocal, pack_bond_info_functor); - copymode = 0; - - k_alocal.modify(); - k_alocal.sync(); -} - -template -KOKKOS_INLINE_FUNCTION -void PairReaxFFKokkos::pack_bond_info_item(const int i) const -{ - const int numbonds = d_numneigh_bonds[i]; - d_alocal(i,0) = tag[i]; - d_alocal(i,1) = type[i]; - d_alocal(i,2) = numbonds; - - int j = 3; - - for (int k = 0; k < numbonds; k++) { - d_alocal(i,j++) = d_neighid(i,k); - } - - d_alocal(i,j++) = molecule.data() ? molecule[i] : 0.0; - - for (int k = 0; k < numbonds; k++) { - d_alocal(i,j++) = d_abo(i,k); - } - - d_alocal(i,j++) = d_total_bo[i]; - d_alocal(i,j++) = paramssing(type[i]).nlp_opt - d_Delta_lp[i]; - d_alocal(i,j++) = q[i]; - - for(; j < d_alocal.extent(1); ++j) { - d_alocal(i,j) = 0.0; - } -} - -/* ---------------------------------------------------------------------- */ - template void PairReaxFFKokkos::FindBondSpecies() { diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index 32007e9970..1ad0955a1e 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -135,7 +135,6 @@ class PairReaxFFKokkos : public PairReaxFF { double memory_usage(); void FindBond(int &); void PackBondBuffer(DAT::tdual_ffloat_1d, int &); - void PackBondInfo(DAT::tdual_float_2d); void FindBondSpecies(); template @@ -293,9 +292,6 @@ class PairReaxFFKokkos : public PairReaxFF { KOKKOS_INLINE_FUNCTION void pack_bond_buffer_item(int, int&, const bool&) const; - KOKKOS_INLINE_FUNCTION - void pack_bond_info_item(const int) const; - KOKKOS_INLINE_FUNCTION void operator()(TagPairReaxFindBondSpeciesZero, const int&) const; @@ -510,8 +506,6 @@ class PairReaxFFKokkos : public PairReaxFF { typename AT::t_ffloat_1d d_buf; DAT::tdual_int_scalar k_nbuf_local; - typename AT::t_float_2d d_alocal; - typedef Kokkos::View t_reax_int4_2d; t_reax_int4_2d d_angular_pack, d_torsion_pack; @@ -555,19 +549,6 @@ struct PairReaxKokkosPackBondBufferFunctor { } }; -template -struct PairReaxKokkosPackBondInfoFunctor { - using device_type = DeviceType; - using value_type = int; - PairReaxFFKokkos c; - PairReaxKokkosPackBondInfoFunctor(PairReaxFFKokkos* c_ptr):c(*c_ptr) {}; - - KOKKOS_INLINE_FUNCTION - void operator()(const int i) const { - c.pack_bond_info_item(i); - } -}; - } #endif From ca143e6ba8da5fba4e44e06472a84ec10fd0f4fe Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Tue, 17 Oct 2023 16:40:03 -0600 Subject: [PATCH 47/60] undo minor change --- examples/reaxff/in.reaxff.tatb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/reaxff/in.reaxff.tatb b/examples/reaxff/in.reaxff.tatb index fd0d846154..f422943a13 100644 --- a/examples/reaxff/in.reaxff.tatb +++ b/examples/reaxff/in.reaxff.tatb @@ -35,7 +35,7 @@ compute bonds all reaxff/bonds variable nqeq equal f_2 # dumps out the local bond information -dump 1 all local 5 bonds_compute.reaxff c_bonds[1] c_bonds[2] c_bonds[3] +dump 1 all local 5 bonds_local.reaxff c_bonds[1] c_bonds[2] c_bonds[3] # dumps out the peratom bond information dump 2 all custom 5 bonds_atom.reaxff c_bonds[*] @@ -57,6 +57,6 @@ timestep 0.0625 # axes yes 0.8 0.02 view 60 -30 #dump_modify 3 pad 3 -#fix 3 all reaxff/species 1 5 5 species.tatb +fix 3 all reaxff/species 1 5 5 species.tatb run 25 From 717e7b064965ec575eeb25ecbbd2d510d9e22100 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Wed, 18 Oct 2023 11:30:16 -0600 Subject: [PATCH 48/60] Address comments --- src/KOKKOS/compute_reaxff_bonds_kokkos.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp b/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp index 921acb9193..5c8885ce7a 100644 --- a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp @@ -39,6 +39,7 @@ ComputeReaxFFBondsKokkos::ComputeReaxFFBondsKokkos(LAMMPS *lmp, int ComputeReaxFFBonds(lmp, narg, arg), nbuf(-1), buf(nullptr) { + kokkosable = 1; } @@ -113,9 +114,8 @@ void ComputeReaxFFBondsKokkos::compute_local() { invoked_local = update->ntimestep; - if(invoked_bonds < update->ntimestep) { + if(invoked_bonds < update->ntimestep) compute_bonds(); - } if(nbonds > prev_nbonds) { // grow array_local @@ -151,9 +151,8 @@ void ComputeReaxFFBondsKokkos::compute_peratom() { invoked_peratom = update->ntimestep; - if(invoked_bonds < update->ntimestep) { + if(invoked_bonds < update->ntimestep) compute_bonds(); - } // extract peratom bond information from buffer int j = 1; From 9bffeb9512b595dd7fd2db4b27c7c46d3e45a034 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Wed, 18 Oct 2023 16:21:10 -0600 Subject: [PATCH 49/60] Next iteration --- examples/reaxff/in.reaxff.tatb | 2 +- src/KOKKOS/compute_reaxff_bonds_kokkos.cpp | 47 +++++++++--------- src/KOKKOS/compute_reaxff_bonds_kokkos.h | 4 +- src/KOKKOS/pair_reaxff_kokkos.cpp | 57 ++++++++++++++++++++++ src/KOKKOS/pair_reaxff_kokkos.h | 17 +++++++ src/REAXFF/compute_reaxff_bonds.cpp | 30 +++++------- src/REAXFF/compute_reaxff_bonds.h | 2 +- 7 files changed, 114 insertions(+), 45 deletions(-) diff --git a/examples/reaxff/in.reaxff.tatb b/examples/reaxff/in.reaxff.tatb index f422943a13..50f572a994 100644 --- a/examples/reaxff/in.reaxff.tatb +++ b/examples/reaxff/in.reaxff.tatb @@ -38,7 +38,7 @@ variable nqeq equal f_2 dump 1 all local 5 bonds_local.reaxff c_bonds[1] c_bonds[2] c_bonds[3] # dumps out the peratom bond information -dump 2 all custom 5 bonds_atom.reaxff c_bonds[*] +dump 2 all custom 5 bonds_atom.reaxff id type q c_bonds[*] thermo 5 thermo_style custom step temp epair etotal press & diff --git a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp b/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp index 5c8885ce7a..f36832e6ac 100644 --- a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp @@ -68,7 +68,7 @@ void ComputeReaxFFBondsKokkos::compute_bonds() if (atom->nlocal > nlocal) { memory->destroy(array_atom); nlocal = atom->nlocal; - memory->create(array_atom, nlocal, 7, "reaxff/bonds:array_atom"); + memory->create(array_atom, nlocal, 3, "reaxff/bonds:array_atom"); } // retrieve bond information from kokkos pair style. the data potentially @@ -76,13 +76,14 @@ void ComputeReaxFFBondsKokkos::compute_bonds() // compute_local and compute_atom then expand the data from this buffer into // appropiate arrays for consumption by others (e.g. dump local, dump custom // or library interface) + int maxnumbonds = 0; if (reaxff->execution_space == Device) device_pair()->FindBond(maxnumbonds); else host_pair()->FindBond(maxnumbonds); - nbuf = 1+(maxnumbonds*2 + 7)*nlocal; + nbuf = (maxnumbonds*2 + 3)*nlocal; if(!buf || k_buf.extent(0) < nbuf) { memoryKK->destroy_kokkos(k_buf, buf); @@ -90,20 +91,21 @@ void ComputeReaxFFBondsKokkos::compute_bonds() } // Pass information to buffer, will sync to host + int nbuf_local; if (reaxff->execution_space == Device) - device_pair()->PackBondBuffer(k_buf, nbuf_local); + device_pair()->PackReducedBondBuffer(k_buf, nbuf_local); else - host_pair()->PackBondBuffer(k_buf, nbuf_local); - buf[0] = nlocal; + host_pair()->PackReducedBondBuffer(k_buf, nbuf_local); // Extract number of bonds from buffer + nbonds = 0; - int j = 1; + int j = 0; for (int i = 0; i < nlocal; i++) { - int numbonds = static_cast(buf[j+5]); + int numbonds = static_cast(buf[j+2]); nbonds += numbonds; - j += 2*numbonds + 7; + j += 2*numbonds + 3; } } @@ -127,20 +129,21 @@ void ComputeReaxFFBondsKokkos::compute_local() size_local_rows = nbonds; // extract local bond information from buffer + int b = 0; - int j = 1; + int j = 0; for (int i = 0; i < nlocal; ++i) { - const int numbonds = static_cast(buf[j+5]); - const int neigh_offset = j + 6; - const int bo_offset = neigh_offset + numbonds + 1; + const int numbonds = static_cast(buf[j+2]); + const int neigh_offset = j + 3; + const int bo_offset = neigh_offset + numbonds; for (int k = 0; k < numbonds; k++) { auto bond = array_local[b++]; bond[0] = i; bond[1] = static_cast (buf[neigh_offset+k]); bond[2] = buf[bo_offset+k]; } - j += 2*numbonds + 7; + j += 2*numbonds + 3; } } @@ -155,19 +158,15 @@ void ComputeReaxFFBondsKokkos::compute_peratom() compute_bonds(); // extract peratom bond information from buffer - int j = 1; + + int j = 0; for (int i = 0; i < nlocal; ++i) { auto ptr = array_atom[i]; - int numbonds = static_cast(buf[j+5]); - const int mol_offset = j + 6 + numbonds; - ptr[0] = buf[j]; // jtag - ptr[1] = buf[j+1]; // itype + int numbonds = static_cast(buf[j+2]); + ptr[0] = buf[j]; // sbo + ptr[1] = buf[j+1]; // nlp ptr[2] = numbonds; - ptr[3] = buf[mol_offset]; // mol - ptr[4] = buf[j+2]; // sbo - ptr[5] = buf[j+3]; // nlp - ptr[6] = buf[j+4]; // q - j += 2*numbonds + 7; + j += 2*numbonds + 3; } } @@ -179,7 +178,7 @@ template double ComputeReaxFFBondsKokkos::memory_usage() { double bytes = (double)(nbonds*3) * sizeof(double); - bytes += (double)(nlocal*7) * sizeof(double); + bytes += (double)(nlocal*3) * sizeof(double); return bytes; } diff --git a/src/KOKKOS/compute_reaxff_bonds_kokkos.h b/src/KOKKOS/compute_reaxff_bonds_kokkos.h index 45020ffa81..48f3860283 100644 --- a/src/KOKKOS/compute_reaxff_bonds_kokkos.h +++ b/src/KOKKOS/compute_reaxff_bonds_kokkos.h @@ -52,11 +52,11 @@ class ComputeReaxFFBondsKokkos : public ComputeReaxFFBonds { typename AT::tdual_float_1d k_buf; auto device_pair() { - return dynamic_cast*>(reaxff); + return static_cast*>(reaxff); } auto host_pair() { - return dynamic_cast*>(reaxff); + return static_cast*>(reaxff); } }; diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index c7d54b80cd..e298eca2da 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4247,6 +4247,30 @@ void PairReaxFFKokkos::PackBondBuffer(DAT::tdual_ffloat_1d k_buf, in nbuf_local = k_nbuf_local.h_view(); } +/* ---------------------------------------------------------------------- */ + +template +void PairReaxFFKokkos::PackReducedBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local) +{ + d_buf = k_buf.view(); + k_params_sing.template sync(); + + copymode = 1; + nlocal = atomKK->nlocal; + PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); + Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + copymode = 0; + + k_buf.modify(); + k_nbuf_local.modify(); + + k_buf.sync(); + k_nbuf_local.sync(); + nbuf_local = k_nbuf_local.h_view(); +} + +/* ---------------------------------------------------------------------- */ + template KOKKOS_INLINE_FUNCTION void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bool &final) const @@ -4288,6 +4312,39 @@ void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bo k_nbuf_local.view()() = j - 1; } +template +KOKKOS_INLINE_FUNCTION +void PairReaxFFKokkos::pack_reduced_bond_buffer_item(int i, int &j, const bool &final) const +{ + const int numbonds = d_numneigh_bonds[i]; + if (final) { + d_buf[j] = d_total_bo[i]; + d_buf[j+1] = paramssing(type[i]).nlp_opt - d_Delta_lp[i]; + d_buf[j+2] = numbonds; + } + + j += 3; + + if (final) { + for (int k = 0; k < numbonds; ++k) { + d_buf[j+k] = d_neighid(i,k); + } + } + + j += numbonds; + + if (final) { + for (int k = 0; k < numbonds; k++) { + d_buf[j+k] = d_abo(i,k); + } + } + + j += numbonds; + + if (final && i == nlocal-1) + k_nbuf_local.view()() = j - 1; +} + /* ---------------------------------------------------------------------- */ template diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index 1ad0955a1e..571dd63fd1 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -135,6 +135,7 @@ class PairReaxFFKokkos : public PairReaxFF { double memory_usage(); void FindBond(int &); void PackBondBuffer(DAT::tdual_ffloat_1d, int &); + void PackReducedBondBuffer(DAT::tdual_ffloat_1d, int &); void FindBondSpecies(); template @@ -292,6 +293,9 @@ class PairReaxFFKokkos : public PairReaxFF { KOKKOS_INLINE_FUNCTION void pack_bond_buffer_item(int, int&, const bool&) const; + KOKKOS_INLINE_FUNCTION + void pack_reduced_bond_buffer_item(int, int&, const bool&) const; + KOKKOS_INLINE_FUNCTION void operator()(TagPairReaxFindBondSpeciesZero, const int&) const; @@ -549,6 +553,19 @@ struct PairReaxKokkosPackBondBufferFunctor { } }; +template +struct PairReaxKokkosPackReducedBondBufferFunctor { + typedef DeviceType device_type; + typedef int value_type; + PairReaxFFKokkos c; + PairReaxKokkosPackReducedBondBufferFunctor(PairReaxFFKokkos* c_ptr):c(*c_ptr) {}; + + KOKKOS_INLINE_FUNCTION + void operator()(const int ii, int &j, const bool &final) const { + c.pack_reduced_bond_buffer_item(ii,j,final); + } +}; + } #endif diff --git a/src/REAXFF/compute_reaxff_bonds.cpp b/src/REAXFF/compute_reaxff_bonds.cpp index 7f1605263c..3b3520fda3 100644 --- a/src/REAXFF/compute_reaxff_bonds.cpp +++ b/src/REAXFF/compute_reaxff_bonds.cpp @@ -35,7 +35,7 @@ using namespace ReaxFF; ComputeReaxFFBonds::ComputeReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - abo(nullptr), neighid(nullptr), numneigh(nullptr), reaxff(nullptr) + abo(nullptr), neighid(nullptr), bondcount(nullptr), reaxff(nullptr) { if (atom->tag_consecutive() == 0) error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds"); @@ -48,7 +48,7 @@ ComputeReaxFFBonds::ComputeReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : nlocal = -1; prev_nbonds = -1; - size_peratom_cols = 7; + size_peratom_cols = 3; size_local_rows = 0; size_local_cols = 3; @@ -65,7 +65,7 @@ ComputeReaxFFBonds::~ComputeReaxFFBonds() memory->destroy(array_atom); memory->destroy(abo); memory->destroy(neighid); - memory->destroy(numneigh); + memory->destroy(bondcount); } /* ---------------------------------------------------------------------- */ @@ -110,7 +110,7 @@ int ComputeReaxFFBonds::FindBond() nj++; } } - numneigh[i] = nj; + bondcount[i] = nj; numbonds += nj; } return numbonds; @@ -126,17 +126,17 @@ void ComputeReaxFFBonds::compute_bonds() if (atom->nlocal > nlocal) { memory->destroy(abo); memory->destroy(neighid); - memory->destroy(numneigh); + memory->destroy(bondcount); memory->destroy(array_atom); nlocal = atom->nlocal; memory->create(abo, nlocal, MAXREAXBOND, "reaxff/bonds:abo"); memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/bonds:neighid"); - memory->create(numneigh, nlocal, "reaxff/bonds:numneigh"); - memory->create(array_atom, nlocal, 7, "reaxff/bonds:array_atom"); + memory->create(bondcount, nlocal, "reaxff/bonds:bondcount"); + memory->create(array_atom, nlocal, 3, "reaxff/bonds:array_atom"); } for (int i = 0; i < nlocal; i++) { - numneigh[i] = 0; + bondcount[i] = 0; for (int j = 0; j < MAXREAXBOND; j++) { neighid[i][j] = 0; abo[i][j] = 0.0; @@ -168,7 +168,7 @@ void ComputeReaxFFBonds::compute_local() int b = 0; for (int i = 0; i < nlocal; ++i) { - const int numbonds = numneigh[i]; + const int numbonds = bondcount[i]; for (int k = 0; k < numbonds; k++) { auto bond = array_local[b++]; @@ -191,13 +191,9 @@ void ComputeReaxFFBonds::compute_peratom() for (int i = 0; i < nlocal; ++i) { auto ptr = array_atom[i]; - ptr[0] = atom->tag[i]; - ptr[1] = atom->type[i]; - ptr[2] = numneigh[i]; - ptr[3] = (atom->molecule == nullptr) ? 0.0 : atom->molecule[i]; - ptr[4] = reaxff->api->workspace->total_bond_order[i]; - ptr[5] = reaxff->api->workspace->nlp[i]; - ptr[6] = atom->q[i]; + ptr[0] = reaxff->api->workspace->total_bond_order[i]; + ptr[1] = reaxff->api->workspace->nlp[i]; + ptr[2] = bondcount[i]; } } @@ -208,7 +204,7 @@ void ComputeReaxFFBonds::compute_peratom() double ComputeReaxFFBonds::memory_usage() { double bytes = (double)(nbonds*3) * sizeof(double); - bytes += (double)(nlocal*7) * sizeof(double); + bytes += (double)(nlocal*3) * sizeof(double); bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); bytes += (double)(nlocal) * sizeof(int); return bytes; diff --git a/src/REAXFF/compute_reaxff_bonds.h b/src/REAXFF/compute_reaxff_bonds.h index 6b00cef7ed..b876c9e02d 100644 --- a/src/REAXFF/compute_reaxff_bonds.h +++ b/src/REAXFF/compute_reaxff_bonds.h @@ -46,7 +46,7 @@ class ComputeReaxFFBonds : public Compute { tagint **neighid; double **abo; - int *numneigh; + int *bondcount; class PairReaxFF *reaxff; private: From afd0107f01ff552aba336ccc02af240f3aed33a8 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Wed, 8 Nov 2023 14:17:18 -0700 Subject: [PATCH 50/60] Add new files to makefile build system --- src/KOKKOS/Install.sh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 489efc55a0..01db058d5b 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -165,6 +165,8 @@ action fix_qeq_reaxff_kokkos.cpp fix_qeq_reaxff.cpp action fix_qeq_reaxff_kokkos.h fix_qeq_reaxff.h action fix_reaxff_bonds_kokkos.cpp fix_reaxff_bonds.cpp action fix_reaxff_bonds_kokkos.h fix_reaxff_bonds.h +action compute_reaxff_bonds_kokkos.cpp compute_reaxff_bonds.cpp +action compute_reaxff_bonds_kokkos.h compute_reaxff_bonds.h action fix_reaxff_species_kokkos.cpp fix_reaxff_species.cpp action fix_reaxff_species_kokkos.h fix_reaxff_species.h action fix_rx_kokkos.cpp fix_rx.cpp From 16f0806da07289e3aea388b453e03ba24f042120 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 20 Nov 2023 15:36:46 -0700 Subject: [PATCH 51/60] Rename compute to reaxff/atom --- examples/reaxff/in.reaxff.tatb | 2 +- src/KOKKOS/Install.sh | 4 +- ...kos.cpp => compute_reaxff_atom_kokkos.cpp} | 46 ++++++------- ..._kokkos.h => compute_reaxff_atom_kokkos.h} | 14 ++-- src/KOKKOS/pair_reaxff_kokkos.cpp | 35 ++++++---- src/KOKKOS/pair_reaxff_kokkos.h | 7 +- ...axff_bonds.cpp => compute_reaxff_atom.cpp} | 64 ++++++++++++------- ...e_reaxff_bonds.h => compute_reaxff_atom.h} | 13 ++-- 8 files changed, 108 insertions(+), 77 deletions(-) rename src/KOKKOS/{compute_reaxff_bonds_kokkos.cpp => compute_reaxff_atom_kokkos.cpp} (76%) rename src/KOKKOS/{compute_reaxff_bonds_kokkos.h => compute_reaxff_atom_kokkos.h} (79%) rename src/REAXFF/{compute_reaxff_bonds.cpp => compute_reaxff_atom.cpp} (74%) rename src/REAXFF/{compute_reaxff_bonds.h => compute_reaxff_atom.h} (84%) diff --git a/examples/reaxff/in.reaxff.tatb b/examples/reaxff/in.reaxff.tatb index 50f572a994..967ed0a1d6 100644 --- a/examples/reaxff/in.reaxff.tatb +++ b/examples/reaxff/in.reaxff.tatb @@ -31,7 +31,7 @@ neigh_modify delay 0 every 5 check no fix 1 all nve fix 2 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff fix 4 all reaxff/bonds 5 bonds.reaxff -compute bonds all reaxff/bonds +compute bonds all reaxff/atom bonds yes variable nqeq equal f_2 # dumps out the local bond information diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 01db058d5b..28fe6d5cd6 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -165,8 +165,8 @@ action fix_qeq_reaxff_kokkos.cpp fix_qeq_reaxff.cpp action fix_qeq_reaxff_kokkos.h fix_qeq_reaxff.h action fix_reaxff_bonds_kokkos.cpp fix_reaxff_bonds.cpp action fix_reaxff_bonds_kokkos.h fix_reaxff_bonds.h -action compute_reaxff_bonds_kokkos.cpp compute_reaxff_bonds.cpp -action compute_reaxff_bonds_kokkos.h compute_reaxff_bonds.h +action compute_reaxff_atom_kokkos.cpp compute_reaxff_atom.cpp +action compute_reaxff_atom_kokkos.h compute_reaxff_atom.h action fix_reaxff_species_kokkos.cpp fix_reaxff_species.cpp action fix_reaxff_species_kokkos.h fix_reaxff_species.h action fix_rx_kokkos.cpp fix_rx.cpp diff --git a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp similarity index 76% rename from src/KOKKOS/compute_reaxff_bonds_kokkos.cpp rename to src/KOKKOS/compute_reaxff_atom_kokkos.cpp index f36832e6ac..2485d9ae72 100644 --- a/src/KOKKOS/compute_reaxff_bonds_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -16,7 +16,7 @@ Contributing author: Richard Berger (LANL) ------------------------------------------------------------------------- */ -#include "compute_reaxff_bonds_kokkos.h" +#include "compute_reaxff_atom_kokkos.h" #include "atom.h" #include "molecule.h" #include "update.h" @@ -35,8 +35,8 @@ using namespace ReaxFF; /* ---------------------------------------------------------------------- */ template -ComputeReaxFFBondsKokkos::ComputeReaxFFBondsKokkos(LAMMPS *lmp, int narg, char **arg) : - ComputeReaxFFBonds(lmp, narg, arg), +ComputeReaxFFAtomKokkos::ComputeReaxFFAtomKokkos(LAMMPS *lmp, int narg, char **arg) : + ComputeReaxFFAtom(lmp, narg, arg), nbuf(-1), buf(nullptr) { kokkosable = 1; @@ -46,7 +46,7 @@ ComputeReaxFFBondsKokkos::ComputeReaxFFBondsKokkos(LAMMPS *lmp, int /* ---------------------------------------------------------------------- */ template -ComputeReaxFFBondsKokkos::~ComputeReaxFFBondsKokkos() +ComputeReaxFFAtomKokkos::~ComputeReaxFFAtomKokkos() { memoryKK->destroy_kokkos(k_buf, buf); } @@ -54,21 +54,21 @@ ComputeReaxFFBondsKokkos::~ComputeReaxFFBondsKokkos() /* ---------------------------------------------------------------------- */ template -void ComputeReaxFFBondsKokkos::init() +void ComputeReaxFFAtomKokkos::init() { reaxff = dynamic_cast(force->pair_match("^reax../kk",0)); - if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/bonds without " + if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/atom without " "pair_style reaxff/kk"); } /* ---------------------------------------------------------------------- */ template -void ComputeReaxFFBondsKokkos::compute_bonds() +void ComputeReaxFFAtomKokkos::compute_bonds() { if (atom->nlocal > nlocal) { memory->destroy(array_atom); nlocal = atom->nlocal; - memory->create(array_atom, nlocal, 3, "reaxff/bonds:array_atom"); + memory->create(array_atom, nlocal, 3, "reaxff/atom:array_atom"); } // retrieve bond information from kokkos pair style. the data potentially @@ -83,20 +83,20 @@ void ComputeReaxFFBondsKokkos::compute_bonds() else host_pair()->FindBond(maxnumbonds); - nbuf = (maxnumbonds*2 + 3)*nlocal; + nbuf = ((store_bonds ? maxnumbonds*2 : 0) + 3)*nlocal; if(!buf || k_buf.extent(0) < nbuf) { memoryKK->destroy_kokkos(k_buf, buf); - memoryKK->create_kokkos(k_buf, buf, nbuf, "reaxff/bonds:buf"); + memoryKK->create_kokkos(k_buf, buf, nbuf, "reaxff/atom:buf"); } // Pass information to buffer, will sync to host int nbuf_local; if (reaxff->execution_space == Device) - device_pair()->PackReducedBondBuffer(k_buf, nbuf_local); + device_pair()->PackReducedBondBuffer(k_buf, nbuf_local, store_bonds); else - host_pair()->PackReducedBondBuffer(k_buf, nbuf_local); + host_pair()->PackReducedBondBuffer(k_buf, nbuf_local, store_bonds); // Extract number of bonds from buffer @@ -105,14 +105,14 @@ void ComputeReaxFFBondsKokkos::compute_bonds() for (int i = 0; i < nlocal; i++) { int numbonds = static_cast(buf[j+2]); nbonds += numbonds; - j += 2*numbonds + 3; + j += (store_bonds ? 2*numbonds : 0) + 3; } } /* ---------------------------------------------------------------------- */ template -void ComputeReaxFFBondsKokkos::compute_local() +void ComputeReaxFFAtomKokkos::compute_local() { invoked_local = update->ntimestep; @@ -122,7 +122,7 @@ void ComputeReaxFFBondsKokkos::compute_local() if(nbonds > prev_nbonds) { // grow array_local memory->destroy(array_local); - memory->create(array_local, nbonds, 3, "reaxff/bonds:array_local"); + memory->create(array_local, nbonds, 3, "reaxff/atom:array_local"); prev_nbonds = nbonds; } @@ -150,7 +150,7 @@ void ComputeReaxFFBondsKokkos::compute_local() /* ---------------------------------------------------------------------- */ template -void ComputeReaxFFBondsKokkos::compute_peratom() +void ComputeReaxFFAtomKokkos::compute_peratom() { invoked_peratom = update->ntimestep; @@ -166,7 +166,7 @@ void ComputeReaxFFBondsKokkos::compute_peratom() ptr[0] = buf[j]; // sbo ptr[1] = buf[j+1]; // nlp ptr[2] = numbonds; - j += 2*numbonds + 3; + j += (store_bonds ? 2*numbonds : 0) + 3; } } @@ -175,16 +175,18 @@ void ComputeReaxFFBondsKokkos::compute_peratom() ------------------------------------------------------------------------- */ template -double ComputeReaxFFBondsKokkos::memory_usage() +double ComputeReaxFFAtomKokkos::memory_usage() { - double bytes = (double)(nbonds*3) * sizeof(double); - bytes += (double)(nlocal*3) * sizeof(double); + double bytes = (double)(nlocal*3) * sizeof(double); + if(store_bonds) + bytes += (double)(nbonds*3) * sizeof(double); + bytes += (double)(nbuf > 0 ? nbuf * sizeof(double) : 0); return bytes; } namespace LAMMPS_NS { -template class ComputeReaxFFBondsKokkos; +template class ComputeReaxFFAtomKokkos; #ifdef LMP_KOKKOS_GPU -template class ComputeReaxFFBondsKokkos; +template class ComputeReaxFFAtomKokkos; #endif } diff --git a/src/KOKKOS/compute_reaxff_bonds_kokkos.h b/src/KOKKOS/compute_reaxff_atom_kokkos.h similarity index 79% rename from src/KOKKOS/compute_reaxff_bonds_kokkos.h rename to src/KOKKOS/compute_reaxff_atom_kokkos.h index 48f3860283..7037c7e308 100644 --- a/src/KOKKOS/compute_reaxff_bonds_kokkos.h +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.h @@ -17,29 +17,29 @@ #ifdef COMPUTE_CLASS // clang-format off -ComputeStyle(reaxff/bonds/kk,ComputeReaxFFBondsKokkos); -ComputeStyle(reaxff/bonds/kk/device,ComputeReaxFFBondsKokkos); -ComputeStyle(reaxff/bonds/kk/host,ComputeReaxFFBondsKokkos); +ComputeStyle(reaxff/atom/kk,ComputeReaxFFAtomKokkos); +ComputeStyle(reaxff/atom/kk/device,ComputeReaxFFAtomKokkos); +ComputeStyle(reaxff/atom/kk/host,ComputeReaxFFAtomKokkos); // clang-format on #else #ifndef LMP_COMPUTE_REAXFF_BONDS_KOKKOS_H #define LMP_COMPUTE_REAXFF_BONDS_KOKKOS_H -#include "compute_reaxff_bonds.h" +#include "compute_reaxff_atom.h" #include "pair_reaxff_kokkos.h" #include "kokkos_type.h" namespace LAMMPS_NS { template -class ComputeReaxFFBondsKokkos : public ComputeReaxFFBonds { +class ComputeReaxFFAtomKokkos : public ComputeReaxFFAtom { public: using device_type = DeviceType; using AT = ArrayTypes; - ComputeReaxFFBondsKokkos(class LAMMPS *, int, char **); - ~ComputeReaxFFBondsKokkos() override; + ComputeReaxFFAtomKokkos(class LAMMPS *, int, char **); + ~ComputeReaxFFAtomKokkos() override; void init() override; void compute_local() override; void compute_peratom() override; diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index e298eca2da..914962f0e6 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4250,15 +4250,21 @@ void PairReaxFFKokkos::PackBondBuffer(DAT::tdual_ffloat_1d k_buf, in /* ---------------------------------------------------------------------- */ template -void PairReaxFFKokkos::PackReducedBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local) +void PairReaxFFKokkos::PackReducedBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local, bool store_bonds) { d_buf = k_buf.view(); k_params_sing.template sync(); copymode = 1; nlocal = atomKK->nlocal; - PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); - Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + if(store_bonds) { + PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); + Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + } else { + PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); + Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + } + copymode = 0; k_buf.modify(); @@ -4313,6 +4319,7 @@ void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bo } template +template KOKKOS_INLINE_FUNCTION void PairReaxFFKokkos::pack_reduced_bond_buffer_item(int i, int &j, const bool &final) const { @@ -4325,21 +4332,23 @@ void PairReaxFFKokkos::pack_reduced_bond_buffer_item(int i, int &j, j += 3; - if (final) { - for (int k = 0; k < numbonds; ++k) { - d_buf[j+k] = d_neighid(i,k); + if constexpr(STORE_BONDS) { + if (final) { + for (int k = 0; k < numbonds; ++k) { + d_buf[j+k] = d_neighid(i,k); + } } - } - j += numbonds; + j += numbonds; - if (final) { - for (int k = 0; k < numbonds; k++) { - d_buf[j+k] = d_abo(i,k); + if (final) { + for (int k = 0; k < numbonds; k++) { + d_buf[j+k] = d_abo(i,k); + } } - } - j += numbonds; + j += numbonds; + } if (final && i == nlocal-1) k_nbuf_local.view()() = j - 1; diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index 571dd63fd1..f246afcc86 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -135,7 +135,7 @@ class PairReaxFFKokkos : public PairReaxFF { double memory_usage(); void FindBond(int &); void PackBondBuffer(DAT::tdual_ffloat_1d, int &); - void PackReducedBondBuffer(DAT::tdual_ffloat_1d, int &); + void PackReducedBondBuffer(DAT::tdual_ffloat_1d, int &, bool); void FindBondSpecies(); template @@ -293,6 +293,7 @@ class PairReaxFFKokkos : public PairReaxFF { KOKKOS_INLINE_FUNCTION void pack_bond_buffer_item(int, int&, const bool&) const; + template KOKKOS_INLINE_FUNCTION void pack_reduced_bond_buffer_item(int, int&, const bool&) const; @@ -553,7 +554,7 @@ struct PairReaxKokkosPackBondBufferFunctor { } }; -template +template struct PairReaxKokkosPackReducedBondBufferFunctor { typedef DeviceType device_type; typedef int value_type; @@ -562,7 +563,7 @@ struct PairReaxKokkosPackReducedBondBufferFunctor { KOKKOS_INLINE_FUNCTION void operator()(const int ii, int &j, const bool &final) const { - c.pack_reduced_bond_buffer_item(ii,j,final); + c.template pack_reduced_bond_buffer_item(ii,j,final); } }; diff --git a/src/REAXFF/compute_reaxff_bonds.cpp b/src/REAXFF/compute_reaxff_atom.cpp similarity index 74% rename from src/REAXFF/compute_reaxff_bonds.cpp rename to src/REAXFF/compute_reaxff_atom.cpp index 3b3520fda3..d975fe42f8 100644 --- a/src/REAXFF/compute_reaxff_bonds.cpp +++ b/src/REAXFF/compute_reaxff_atom.cpp @@ -16,7 +16,7 @@ Contributing author: Richard Berger (LANL) ------------------------------------------------------------------------- */ -#include "compute_reaxff_bonds.h" +#include "compute_reaxff_atom.h" #include "atom.h" #include "molecule.h" #include "update.h" @@ -33,14 +33,13 @@ using namespace ReaxFF; /* ---------------------------------------------------------------------- */ -ComputeReaxFFBonds::ComputeReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : +ComputeReaxFFAtom::ComputeReaxFFAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), abo(nullptr), neighid(nullptr), bondcount(nullptr), reaxff(nullptr) { if (atom->tag_consecutive() == 0) - error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/bonds"); + error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/atom"); - local_flag = 1; peratom_flag = 1; // initialize output @@ -54,12 +53,25 @@ ComputeReaxFFBonds::ComputeReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : size_local_cols = 3; invoked_bonds = -1; + + store_bonds = false; + + int iarg = 3; + while (iarg narg) utils::missing_cmd_args(FLERR, "compute reaxff/atom bonds", error); + store_bonds = utils::logical(FLERR, arg[iarg+1], false, lmp); + iarg += 2; + } else error->all(FLERR,"Illegal compute reaxff/atom command"); + } + + local_flag = store_bonds; } /* ---------------------------------------------------------------------- */ -ComputeReaxFFBonds::~ComputeReaxFFBonds() +ComputeReaxFFAtom::~ComputeReaxFFAtom() { memory->destroy(array_local); memory->destroy(array_atom); @@ -70,16 +82,16 @@ ComputeReaxFFBonds::~ComputeReaxFFBonds() /* ---------------------------------------------------------------------- */ -void ComputeReaxFFBonds::init() +void ComputeReaxFFAtom::init() { reaxff = dynamic_cast(force->pair_match("^reax..",0)); - if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/bonds without " + if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/atom without " "pair_style reaxff, reaxff/kk, or reaxff/omp"); } /* ---------------------------------------------------------------------- */ -int ComputeReaxFFBonds::FindBond() +int ComputeReaxFFAtom::FindBond() { int *ilist, i, ii, inum; int j, pj, nj; @@ -105,8 +117,10 @@ int ComputeReaxFFBonds::FindBond() bo_tmp = bo_ij->bo_data.BO; if (bo_tmp > bo_cut) { - neighid[i][nj] = jtag; - abo[i][nj] = bo_tmp; + if(store_bonds) { + neighid[i][nj] = jtag; + abo[i][nj] = bo_tmp; + } nj++; } } @@ -119,7 +133,7 @@ int ComputeReaxFFBonds::FindBond() /* ---------------------------------------------------------------------- */ -void ComputeReaxFFBonds::compute_bonds() +void ComputeReaxFFAtom::compute_bonds() { invoked_bonds = update->ntimestep; @@ -129,15 +143,17 @@ void ComputeReaxFFBonds::compute_bonds() memory->destroy(bondcount); memory->destroy(array_atom); nlocal = atom->nlocal; - memory->create(abo, nlocal, MAXREAXBOND, "reaxff/bonds:abo"); - memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/bonds:neighid"); - memory->create(bondcount, nlocal, "reaxff/bonds:bondcount"); - memory->create(array_atom, nlocal, 3, "reaxff/bonds:array_atom"); + if(store_bonds) { + memory->create(abo, nlocal, MAXREAXBOND, "reaxff/atom:abo"); + memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/atom:neighid"); + } + memory->create(bondcount, nlocal, "reaxff/atom:bondcount"); + memory->create(array_atom, nlocal, 3, "reaxff/atom:array_atom"); } for (int i = 0; i < nlocal; i++) { bondcount[i] = 0; - for (int j = 0; j < MAXREAXBOND; j++) { + for (int j = 0; store_bonds && j < MAXREAXBOND; j++) { neighid[i][j] = 0; abo[i][j] = 0.0; } @@ -148,7 +164,7 @@ void ComputeReaxFFBonds::compute_bonds() /* ---------------------------------------------------------------------- */ -void ComputeReaxFFBonds::compute_local() +void ComputeReaxFFAtom::compute_local() { invoked_local = update->ntimestep; @@ -159,7 +175,7 @@ void ComputeReaxFFBonds::compute_local() if(nbonds > prev_nbonds) { // grow array_local memory->destroy(array_local); - memory->create(array_local, nbonds, 3, "reaxff/bonds:array_local"); + memory->create(array_local, nbonds, 3, "reaxff/atom:array_local"); prev_nbonds = nbonds; } @@ -181,7 +197,7 @@ void ComputeReaxFFBonds::compute_local() /* ---------------------------------------------------------------------- */ -void ComputeReaxFFBonds::compute_peratom() +void ComputeReaxFFAtom::compute_peratom() { invoked_peratom = update->ntimestep; @@ -201,11 +217,13 @@ void ComputeReaxFFBonds::compute_peratom() memory usage of local data ------------------------------------------------------------------------- */ -double ComputeReaxFFBonds::memory_usage() +double ComputeReaxFFAtom::memory_usage() { - double bytes = (double)(nbonds*3) * sizeof(double); - bytes += (double)(nlocal*3) * sizeof(double); - bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); + double bytes = (double)(nlocal*3) * sizeof(double); bytes += (double)(nlocal) * sizeof(int); + if(store_bonds) { + bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); + bytes += (double)(nbonds*3) * sizeof(double); + } return bytes; } diff --git a/src/REAXFF/compute_reaxff_bonds.h b/src/REAXFF/compute_reaxff_atom.h similarity index 84% rename from src/REAXFF/compute_reaxff_bonds.h rename to src/REAXFF/compute_reaxff_atom.h index b876c9e02d..31b18e7238 100644 --- a/src/REAXFF/compute_reaxff_bonds.h +++ b/src/REAXFF/compute_reaxff_atom.h @@ -17,21 +17,21 @@ #ifdef COMPUTE_CLASS // clang-format off -ComputeStyle(reaxff/bonds,ComputeReaxFFBonds); +ComputeStyle(reaxff/atom,ComputeReaxFFAtom); // clang-format on #else -#ifndef LMP_COMPUTE_REAXFF_BONDS_H -#define LMP_COMPUTE_REAXFF_BONDS_H +#ifndef LMP_COMPUTE_REAXFF_ATOM_H +#define LMP_COMPUTE_REAXFF_ATOM_H #include "compute.h" namespace LAMMPS_NS { -class ComputeReaxFFBonds : public Compute { +class ComputeReaxFFAtom : public Compute { public: - ComputeReaxFFBonds(class LAMMPS *, int, char **); - ~ComputeReaxFFBonds() override; + ComputeReaxFFAtom(class LAMMPS *, int, char **); + ~ComputeReaxFFAtom() override; void init() override; void compute_local() override; void compute_peratom() override; @@ -43,6 +43,7 @@ class ComputeReaxFFBonds : public Compute { int nlocal; int nbonds; int prev_nbonds; + bool store_bonds; tagint **neighid; double **abo; From fd83ed4004a8d26b95f1f69e75b288804361bfd5 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Tue, 21 Nov 2023 13:42:51 -0700 Subject: [PATCH 52/60] compute reaxff/atom: add support for pair hybrid --- src/KOKKOS/compute_reaxff_atom_kokkos.cpp | 9 ++++--- src/REAXFF/compute_reaxff_atom.cpp | 30 ++++++++++++++++++++--- src/REAXFF/compute_reaxff_atom.h | 1 + 3 files changed, 33 insertions(+), 7 deletions(-) diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp index 2485d9ae72..845e8bc6e7 100644 --- a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -56,9 +56,12 @@ ComputeReaxFFAtomKokkos::~ComputeReaxFFAtomKokkos() template void ComputeReaxFFAtomKokkos::init() { - reaxff = dynamic_cast(force->pair_match("^reax../kk",0)); - if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/atom without " - "pair_style reaxff/kk"); + ComputeReaxFFAtom::init(); + + if(!reaxff || !reaxff->kokkosable) { + error->all(FLERR,"Cannot use compute reaxff/atom/kk without " + "pair_style reaxff/kk"); + } } /* ---------------------------------------------------------------------- */ diff --git a/src/REAXFF/compute_reaxff_atom.cpp b/src/REAXFF/compute_reaxff_atom.cpp index d975fe42f8..1d7bdf0e7d 100644 --- a/src/REAXFF/compute_reaxff_atom.cpp +++ b/src/REAXFF/compute_reaxff_atom.cpp @@ -55,10 +55,21 @@ ComputeReaxFFAtom::ComputeReaxFFAtom(LAMMPS *lmp, int narg, char **arg) : invoked_bonds = -1; store_bonds = false; + nsub = 0; int iarg = 3; while (iarg narg) utils::missing_cmd_args(FLERR, "compute reaxff/atom pair", error); + ++iarg; + + if (isdigit(arg[iarg][0])) { + nsub = utils::inumeric(FLERR, arg[iarg], false, lmp); + ++iarg; + if (nsub > 0) continue; + } + error->all(FLERR, "Illegal compute reaxff/atom command"); + } else if (strcmp(arg[iarg], "bonds") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "compute reaxff/atom bonds", error); store_bonds = utils::logical(FLERR, arg[iarg+1], false, lmp); iarg += 2; @@ -84,9 +95,20 @@ ComputeReaxFFAtom::~ComputeReaxFFAtom() void ComputeReaxFFAtom::init() { - reaxff = dynamic_cast(force->pair_match("^reax..",0)); - if (reaxff == nullptr) error->all(FLERR,"Cannot use compute reaxff/atom without " - "pair_style reaxff, reaxff/kk, or reaxff/omp"); + if (lmp->suffix_enable) { + if (lmp->suffix) + reaxff = dynamic_cast(force->pair_match(fmt::format("^reax../{}", lmp->suffix), 0, nsub)); + if (!reaxff && lmp->suffix2) + reaxff = dynamic_cast(force->pair_match(fmt::format("^reax../{}", lmp->suffix2), 0, nsub)); + } + + if (!reaxff) reaxff = dynamic_cast(force->pair_match("^reax..", 0, nsub)); + + if (!reaxff) error->all(FLERR,"Cannot use compute reaxff/atom without " + "pair_style reaxff or reaxff/omp"); + + if(reaxff->kokkosable && !kokkosable) + error->all(FLERR,"Cannot use compute reaxff/atom with pair_style reaxff/kk. Use reaxff/atom/kk."); } /* ---------------------------------------------------------------------- */ diff --git a/src/REAXFF/compute_reaxff_atom.h b/src/REAXFF/compute_reaxff_atom.h index 31b18e7238..1f9aaec1ae 100644 --- a/src/REAXFF/compute_reaxff_atom.h +++ b/src/REAXFF/compute_reaxff_atom.h @@ -43,6 +43,7 @@ class ComputeReaxFFAtom : public Compute { int nlocal; int nbonds; int prev_nbonds; + int nsub; bool store_bonds; tagint **neighid; From a5cc181358cab34f2edfa1ab79e04d71723bfa70 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Tue, 21 Nov 2023 14:03:10 -0700 Subject: [PATCH 53/60] Start with compute reaxff/atom documentation --- doc/src/Commands_compute.rst | 1 + doc/src/compute_reaxff_atom.rst | 67 +++++++++++++++++++++++++++++++++ doc/src/pair_reaxff.rst | 3 +- 3 files changed, 70 insertions(+), 1 deletion(-) create mode 100644 doc/src/compute_reaxff_atom.rst diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index dbd6b58ce7..819a1b10d8 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -115,6 +115,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`property/grid ` * :doc:`property/local ` * :doc:`ptm/atom ` + * :doc:`reaxff/atom (k) ` * :doc:`rdf ` * :doc:`reduce ` * :doc:`reduce/chunk ` diff --git a/doc/src/compute_reaxff_atom.rst b/doc/src/compute_reaxff_atom.rst new file mode 100644 index 0000000000..4906f3fe5c --- /dev/null +++ b/doc/src/compute_reaxff_atom.rst @@ -0,0 +1,67 @@ +.. index:: fix reaxff/atom +.. index:: fix reaxff/atom/kk + +fix reaxff/atom command +======================= + +Accelerator Variants: *reaxff/atom/kk* + +Syntax +"""""" + +.. code-block:: LAMMPS + + compute ID group-ID reaxff/atom attribute args ... keyword value ... + +* ID, group-ID are documented in :doc:`compute ` command +* reaxff/atom = name of this compute command +* attribute = *pair* + + .. parsed-literal:: + + *pair* args = nsub + nsub = *n*-instance of a sub-style, if a pair style is used multiple times in a hybrid style + +* keyword = *bonds* + + .. parsed-literal:: + + *bonds* value = *no* or *yes* + *no* = ignore list of local bonds + *yes* = include list of local bonds + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all reaxff/atom bonds yes + +Description +""""""""""" + +TODO + +---------- + +.. include:: accel_styles.rst + +---------- + +Restrictions +"""""""""""" + +The compute reaxff/atom command requires that the :doc:`pair_style reaxff +` is invoked. This fix is part of the REAXFF package. It is only +enabled if LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +Related commands +"""""""""""""""" + +:doc:`pair_style reaxff ` + +Default +""""""" + +The option defaults are *bonds* = *no*. diff --git a/doc/src/pair_reaxff.rst b/doc/src/pair_reaxff.rst index d28e15b0a2..03d53d1ff4 100644 --- a/doc/src/pair_reaxff.rst +++ b/doc/src/pair_reaxff.rst @@ -373,7 +373,8 @@ Related commands :doc:`pair_coeff `, :doc:`fix qeq/reaxff `, :doc:`fix acks2/reaxff `, :doc:`fix reaxff/bonds `, -:doc:`fix reaxff/species ` +:doc:`fix reaxff/species `, +:doc:`compute reaxff/atom ` Default """"""" From b72c34d4979b184adbb4aea83240e7bcb30e040e Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Tue, 21 Nov 2023 14:20:17 -0700 Subject: [PATCH 54/60] compute reaxff/atom: return tag[i] instead of i --- src/KOKKOS/compute_reaxff_atom_kokkos.cpp | 3 ++- src/REAXFF/compute_reaxff_atom.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp index 845e8bc6e7..c703bc3552 100644 --- a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -135,6 +135,7 @@ void ComputeReaxFFAtomKokkos::compute_local() int b = 0; int j = 0; + auto tag = atom->tag; for (int i = 0; i < nlocal; ++i) { const int numbonds = static_cast(buf[j+2]); @@ -142,7 +143,7 @@ void ComputeReaxFFAtomKokkos::compute_local() const int bo_offset = neigh_offset + numbonds; for (int k = 0; k < numbonds; k++) { auto bond = array_local[b++]; - bond[0] = i; + bond[0] = tag[i]; bond[1] = static_cast (buf[neigh_offset+k]); bond[2] = buf[bo_offset+k]; } diff --git a/src/REAXFF/compute_reaxff_atom.cpp b/src/REAXFF/compute_reaxff_atom.cpp index 1d7bdf0e7d..2e5d4058b8 100644 --- a/src/REAXFF/compute_reaxff_atom.cpp +++ b/src/REAXFF/compute_reaxff_atom.cpp @@ -202,6 +202,7 @@ void ComputeReaxFFAtom::compute_local() } size_local_rows = nbonds; + auto tag = atom->tag; int b = 0; @@ -210,7 +211,7 @@ void ComputeReaxFFAtom::compute_local() for (int k = 0; k < numbonds; k++) { auto bond = array_local[b++]; - bond[0] = i; + bond[0] = tag[i]; bond[1] = neighid[i][k]; bond[2] = abo[i][k]; } From e241f08cfe3c24a49d46cf9e7515b4db335406e4 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 11 Dec 2023 14:05:42 -0700 Subject: [PATCH 55/60] Finish first version of compute reaxff/atom docs --- doc/src/compute.rst | 1 + doc/src/compute_reaxff_atom.rst | 38 ++++++++++++++++++++++++++++----- 2 files changed, 34 insertions(+), 5 deletions(-) diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6737203618..6ef093c16d 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -280,6 +280,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`property/local ` - convert local attributes to local vectors/arrays * :doc:`ptm/atom ` - determines the local lattice structure based on the Polyhedral Template Matching method * :doc:`rdf ` - radial distribution function :math:`g(r)` histogram of group of atoms +* :doc:`reaxff/atom ` - extract ReaxFF bond information * :doc:`reduce ` - combine per-atom quantities into a single global value * :doc:`reduce/chunk ` - reduce per-atom quantities within each chunk * :doc:`reduce/region ` - same as compute reduce, within a region diff --git a/doc/src/compute_reaxff_atom.rst b/doc/src/compute_reaxff_atom.rst index 4906f3fe5c..9f690d19c3 100644 --- a/doc/src/compute_reaxff_atom.rst +++ b/doc/src/compute_reaxff_atom.rst @@ -1,8 +1,8 @@ -.. index:: fix reaxff/atom -.. index:: fix reaxff/atom/kk +.. index:: compute reaxff/atom +.. index:: compute reaxff/atom/kk -fix reaxff/atom command -======================= +compute reaxff/atom command +=========================== Accelerator Variants: *reaxff/atom/kk* @@ -40,7 +40,35 @@ Examples Description """"""""""" -TODO +Define a computation that extracts bond information computed by the ReaxFF +potential specified by :doc:`pair_style reaxff `. + +By default, it produces per-atom data that includes the following columns: + +* abo = atom bond order (sum of all bonds) +* nlp = number of lone pairs +* nb = number of bonds + +Bonds will only be included if its atoms are in the group. + +In addition, if ``bonds`` is set to ``yes``, the compute will also produce a +local array of all bonds on the current processor whose atoms are in the group. +The columns of each entry of this local array are: + +* id_i = atom i id of bond +* id_j = atom j id of bond +* bo = bond order of bond + +Output info +""""""""""" + +This compute calculates a per-atom array and local array depending on the +number of keywords. The number of rows in the local array is the number of +bonds as described above. Both per-atom and local array have 3 columns. + +The arrays can be accessed by any command that uses local and per-atom values +from a compute as input. See the :doc:`Howto output ` page for +an overview of LAMMPS output options. ---------- From 930fbe8c5df277cb0aa1195a5a7ad031a2585a5c Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 11 Dec 2023 14:07:48 -0700 Subject: [PATCH 56/60] reaxff/atom: First attempt of filtering by group --- src/KOKKOS/compute_reaxff_atom_kokkos.cpp | 4 +-- src/KOKKOS/pair_reaxff_kokkos.cpp | 37 +++++++++++++---------- src/KOKKOS/pair_reaxff_kokkos.h | 10 +++--- src/REAXFF/compute_reaxff_atom.cpp | 31 +++++++++++-------- 4 files changed, 47 insertions(+), 35 deletions(-) diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp index c703bc3552..c3e9367ff4 100644 --- a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -82,9 +82,9 @@ void ComputeReaxFFAtomKokkos::compute_bonds() int maxnumbonds = 0; if (reaxff->execution_space == Device) - device_pair()->FindBond(maxnumbonds); + device_pair()->FindBond(maxnumbonds, groupbit); else - host_pair()->FindBond(maxnumbonds); + host_pair()->FindBond(maxnumbonds, groupbit); nbuf = ((store_bonds ? maxnumbonds*2 : 0) + 3)*nlocal; diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index 914962f0e6..29ac398128 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4162,22 +4162,23 @@ double PairReaxFFKokkos::memory_usage() /* ---------------------------------------------------------------------- */ template -void PairReaxFFKokkos::FindBond(int &numbonds) +void PairReaxFFKokkos::FindBond(int &numbonds, int groupbit) { copymode = 1; Kokkos::parallel_for(Kokkos::RangePolicy(0,nmax),*this); bo_cut_bond = api->control->bg_cut; - atomKK->sync(execution_space,TAG_MASK); + atomKK->sync(execution_space,TAG_MASK|MASK_MASK); tag = atomKK->k_tag.view(); + mask = atomKK->k_mask.view(); const int inum = list->inum; NeighListKokkos* k_list = static_cast*>(list); d_ilist = k_list->d_ilist; numbonds = 0; - PairReaxKokkosFindBondFunctor find_bond_functor(this); + PairReaxKokkosFindBondFunctor find_bond_functor(this, groupbit); Kokkos::parallel_reduce(inum,find_bond_functor,numbonds); copymode = 0; } @@ -4194,24 +4195,28 @@ void PairReaxFFKokkos::operator()(TagPairReaxFindBondZero, const int template KOKKOS_INLINE_FUNCTION -void PairReaxFFKokkos::calculate_find_bond_item(int ii, int &numbonds) const +void PairReaxFFKokkos::calculate_find_bond_item(int ii, int &numbonds, int groupbit) const { const int i = d_ilist[ii]; int nj = 0; - const int j_start = d_bo_first[i]; - const int j_end = j_start + d_bo_num[i]; - for (int jj = j_start; jj < j_end; jj++) { - int j = d_bo_list[jj]; - j &= NEIGHMASK; - const tagint jtag = tag[j]; - const int j_index = jj - j_start; - double bo_tmp = d_BO(i,j_index); + if(mask[i] & groupbit) { + const int j_start = d_bo_first[i]; + const int j_end = j_start + d_bo_num[i]; + for (int jj = j_start; jj < j_end; jj++) { + int j = d_bo_list[jj]; + j &= NEIGHMASK; + if(mask[j] & groupbit) { + const tagint jtag = tag[j]; + const int j_index = jj - j_start; + double bo_tmp = d_BO(i,j_index); - if (bo_tmp > bo_cut_bond) { - d_neighid(i,nj) = jtag; - d_abo(i,nj) = bo_tmp; - nj++; + if (bo_tmp > bo_cut_bond) { + d_neighid(i,nj) = jtag; + d_abo(i,nj) = bo_tmp; + nj++; + } + } } } d_numneigh_bonds[i] = nj; diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index f246afcc86..f54350b04b 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -133,7 +133,7 @@ class PairReaxFFKokkos : public PairReaxFF { void compute(int, int); void init_style(); double memory_usage(); - void FindBond(int &); + void FindBond(int &, int groupbit = 1); void PackBondBuffer(DAT::tdual_ffloat_1d, int &); void PackReducedBondBuffer(DAT::tdual_ffloat_1d, int &, bool); void FindBondSpecies(); @@ -288,7 +288,7 @@ class PairReaxFFKokkos : public PairReaxFF { void operator()(TagPairReaxFindBondZero, const int&) const; KOKKOS_INLINE_FUNCTION - void calculate_find_bond_item(int, int&) const; + void calculate_find_bond_item(int, int&, int) const; KOKKOS_INLINE_FUNCTION void pack_bond_buffer_item(int, int&, const bool&) const; @@ -417,6 +417,7 @@ class PairReaxFFKokkos : public PairReaxFF { typename AT::t_f_array f; typename AT::t_int_1d_randomread type; typename AT::t_tagint_1d_randomread tag; + typename AT::t_int_1d_randomread mask; typename AT::t_float_1d_randomread q; typename AT::t_tagint_1d_randomread molecule; @@ -526,8 +527,9 @@ template struct PairReaxKokkosFindBondFunctor { typedef DeviceType device_type; typedef int value_type; + int groupbit; PairReaxFFKokkos c; - PairReaxKokkosFindBondFunctor(PairReaxFFKokkos* c_ptr):c(*c_ptr) {}; + PairReaxKokkosFindBondFunctor(PairReaxFFKokkos* c_ptr, int groupbit):c(*c_ptr),groupbit(groupbit) {}; KOKKOS_INLINE_FUNCTION void join(int &dst, @@ -537,7 +539,7 @@ struct PairReaxKokkosFindBondFunctor { KOKKOS_INLINE_FUNCTION void operator()(const int ii, int &numbonds) const { - c.calculate_find_bond_item(ii,numbonds); + c.calculate_find_bond_item(ii,numbonds,groupbit); } }; diff --git a/src/REAXFF/compute_reaxff_atom.cpp b/src/REAXFF/compute_reaxff_atom.cpp index 2e5d4058b8..fba67ecacf 100644 --- a/src/REAXFF/compute_reaxff_atom.cpp +++ b/src/REAXFF/compute_reaxff_atom.cpp @@ -126,28 +126,33 @@ int ComputeReaxFFAtom::FindBond() bo_cut = reaxff->api->control->bg_cut; tagint *tag = atom->tag; + int * mask = atom->mask; int numbonds = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - nj = 0; + if (mask[i] & groupbit) { + nj = 0; - for (pj = Start_Index(i, reaxff->api->lists); pj < End_Index(i, reaxff->api->lists); ++pj) { - bo_ij = &(reaxff->api->lists->select.bond_list[pj]); - j = bo_ij->nbr; - jtag = tag[j]; - bo_tmp = bo_ij->bo_data.BO; + for (pj = Start_Index(i, reaxff->api->lists); pj < End_Index(i, reaxff->api->lists); ++pj) { + bo_ij = &(reaxff->api->lists->select.bond_list[pj]); + j = bo_ij->nbr; + if (mask[j] & groupbit) { + jtag = tag[j]; + bo_tmp = bo_ij->bo_data.BO; - if (bo_tmp > bo_cut) { - if(store_bonds) { - neighid[i][nj] = jtag; - abo[i][nj] = bo_tmp; + if (bo_tmp > bo_cut) { + if(store_bonds) { + neighid[i][nj] = jtag; + abo[i][nj] = bo_tmp; + } + nj++; + } } - nj++; } + bondcount[i] = nj; + numbonds += nj; } - bondcount[i] = nj; - numbonds += nj; } return numbonds; } From 0cf4c9e7a32585ca7c78d4fd4df4f50c89d83d9b Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Wed, 13 Dec 2023 11:07:42 -0700 Subject: [PATCH 57/60] Whitespace --- src/KOKKOS/compute_reaxff_atom_kokkos.cpp | 13 ++++++------- src/KOKKOS/pair_reaxff_kokkos.cpp | 6 +++--- src/REAXFF/compute_reaxff_atom.cpp | 17 +++++++---------- 3 files changed, 16 insertions(+), 20 deletions(-) diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp index c3e9367ff4..8dbcb9441e 100644 --- a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -42,7 +42,6 @@ ComputeReaxFFAtomKokkos::ComputeReaxFFAtomKokkos(LAMMPS *lmp, int na kokkosable = 1; } - /* ---------------------------------------------------------------------- */ template @@ -58,7 +57,7 @@ void ComputeReaxFFAtomKokkos::init() { ComputeReaxFFAtom::init(); - if(!reaxff || !reaxff->kokkosable) { + if (!reaxff || !reaxff->kokkosable) { error->all(FLERR,"Cannot use compute reaxff/atom/kk without " "pair_style reaxff/kk"); } @@ -88,7 +87,7 @@ void ComputeReaxFFAtomKokkos::compute_bonds() nbuf = ((store_bonds ? maxnumbonds*2 : 0) + 3)*nlocal; - if(!buf || k_buf.extent(0) < nbuf) { + if (!buf || k_buf.extent(0) < nbuf) { memoryKK->destroy_kokkos(k_buf, buf); memoryKK->create_kokkos(k_buf, buf, nbuf, "reaxff/atom:buf"); } @@ -119,10 +118,10 @@ void ComputeReaxFFAtomKokkos::compute_local() { invoked_local = update->ntimestep; - if(invoked_bonds < update->ntimestep) + if (invoked_bonds < update->ntimestep) compute_bonds(); - if(nbonds > prev_nbonds) { + if (nbonds > prev_nbonds) { // grow array_local memory->destroy(array_local); memory->create(array_local, nbonds, 3, "reaxff/atom:array_local"); @@ -158,7 +157,7 @@ void ComputeReaxFFAtomKokkos::compute_peratom() { invoked_peratom = update->ntimestep; - if(invoked_bonds < update->ntimestep) + if (invoked_bonds < update->ntimestep) compute_bonds(); // extract peratom bond information from buffer @@ -182,7 +181,7 @@ template double ComputeReaxFFAtomKokkos::memory_usage() { double bytes = (double)(nlocal*3) * sizeof(double); - if(store_bonds) + if (store_bonds) bytes += (double)(nbonds*3) * sizeof(double); bytes += (double)(nbuf > 0 ? nbuf * sizeof(double) : 0); return bytes; diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index 29ac398128..11a40970c2 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4200,13 +4200,13 @@ void PairReaxFFKokkos::calculate_find_bond_item(int ii, int &numbond const int i = d_ilist[ii]; int nj = 0; - if(mask[i] & groupbit) { + if (mask[i] & groupbit) { const int j_start = d_bo_first[i]; const int j_end = j_start + d_bo_num[i]; for (int jj = j_start; jj < j_end; jj++) { int j = d_bo_list[jj]; j &= NEIGHMASK; - if(mask[j] & groupbit) { + if (mask[j] & groupbit) { const tagint jtag = tag[j]; const int j_index = jj - j_start; double bo_tmp = d_BO(i,j_index); @@ -4262,7 +4262,7 @@ void PairReaxFFKokkos::PackReducedBondBuffer(DAT::tdual_ffloat_1d k_ copymode = 1; nlocal = atomKK->nlocal; - if(store_bonds) { + if (store_bonds) { PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); } else { diff --git a/src/REAXFF/compute_reaxff_atom.cpp b/src/REAXFF/compute_reaxff_atom.cpp index fba67ecacf..1b54f3b82a 100644 --- a/src/REAXFF/compute_reaxff_atom.cpp +++ b/src/REAXFF/compute_reaxff_atom.cpp @@ -79,7 +79,6 @@ ComputeReaxFFAtom::ComputeReaxFFAtom(LAMMPS *lmp, int narg, char **arg) : local_flag = store_bonds; } - /* ---------------------------------------------------------------------- */ ComputeReaxFFAtom::~ComputeReaxFFAtom() @@ -107,7 +106,7 @@ void ComputeReaxFFAtom::init() if (!reaxff) error->all(FLERR,"Cannot use compute reaxff/atom without " "pair_style reaxff or reaxff/omp"); - if(reaxff->kokkosable && !kokkosable) + if (reaxff->kokkosable && !kokkosable) error->all(FLERR,"Cannot use compute reaxff/atom with pair_style reaxff/kk. Use reaxff/atom/kk."); } @@ -142,7 +141,7 @@ int ComputeReaxFFAtom::FindBond() bo_tmp = bo_ij->bo_data.BO; if (bo_tmp > bo_cut) { - if(store_bonds) { + if (store_bonds) { neighid[i][nj] = jtag; abo[i][nj] = bo_tmp; } @@ -157,7 +156,6 @@ int ComputeReaxFFAtom::FindBond() return numbonds; } - /* ---------------------------------------------------------------------- */ void ComputeReaxFFAtom::compute_bonds() @@ -170,7 +168,7 @@ void ComputeReaxFFAtom::compute_bonds() memory->destroy(bondcount); memory->destroy(array_atom); nlocal = atom->nlocal; - if(store_bonds) { + if (store_bonds) { memory->create(abo, nlocal, MAXREAXBOND, "reaxff/atom:abo"); memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/atom:neighid"); } @@ -195,11 +193,10 @@ void ComputeReaxFFAtom::compute_local() { invoked_local = update->ntimestep; - if(invoked_bonds < update->ntimestep) { + if (invoked_bonds < update->ntimestep) compute_bonds(); - } - if(nbonds > prev_nbonds) { + if (nbonds > prev_nbonds) { // grow array_local memory->destroy(array_local); memory->create(array_local, nbonds, 3, "reaxff/atom:array_local"); @@ -229,7 +226,7 @@ void ComputeReaxFFAtom::compute_peratom() { invoked_peratom = update->ntimestep; - if(invoked_bonds < update->ntimestep) { + if (invoked_bonds < update->ntimestep) { compute_bonds(); } @@ -249,7 +246,7 @@ double ComputeReaxFFAtom::memory_usage() { double bytes = (double)(nlocal*3) * sizeof(double); bytes += (double)(nlocal) * sizeof(int); - if(store_bonds) { + if (store_bonds) { bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); bytes += (double)(nbonds*3) * sizeof(double); } From 603837c96ca07d84c638b9db1cb1f39b9cddbc0d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 14 Dec 2023 16:21:04 -0500 Subject: [PATCH 58/60] add versionadded tag --- doc/src/compute_reaxff_atom.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/src/compute_reaxff_atom.rst b/doc/src/compute_reaxff_atom.rst index 9f690d19c3..997ad02e9f 100644 --- a/doc/src/compute_reaxff_atom.rst +++ b/doc/src/compute_reaxff_atom.rst @@ -40,6 +40,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Define a computation that extracts bond information computed by the ReaxFF potential specified by :doc:`pair_style reaxff `. From 7dab2b7eee941e28e01416ede9688e661c1cf520 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 14 Dec 2023 16:21:17 -0500 Subject: [PATCH 59/60] add new package files to .gitignore --- src/.gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/.gitignore b/src/.gitignore index 3b7902303d..112a1486f7 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -633,6 +633,8 @@ /compute_ptm_atom.h /compute_rattlers_atom.cpp /compute_rattlers_atom.h +/compute_reaxff_atom.cpp +/compute_reaxff_atom.h /compute_rigid_local.cpp /compute_rigid_local.h /compute_slcsa_atom.cpp From 1df91f21a18d88006e07ec50e56c16ca7036b927 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 14 Dec 2023 17:59:18 -0500 Subject: [PATCH 60/60] workaround hack for macOS --- src/PYTHON/python_impl.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index c1c26c7bb3..ea4ac63ce7 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -79,7 +79,7 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp) // Force the stdout and stderr streams to be unbuffered. bool unbuffered = PYTHONUNBUFFERED != nullptr && strcmp(PYTHONUNBUFFERED, "1") == 0; -#if PY_VERSION_HEX >= 0x030800f0 +#if (PY_VERSION_HEX >= 0x030800f0) && !defined(__APPLE__) PyConfig config; PyConfig_InitPythonConfig(&config); config.buffered_stdio = !unbuffered; @@ -87,6 +87,8 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp) // Python Global configuration variable Py_UnbufferedStdioFlag = unbuffered; #endif +#else +#warning Cannot force stdout and stderr to be unbuffered #endif #ifdef MLIAP_PYTHON