From 541680c79804a67488d48d95a6f7c0feb882e301 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 4 May 2024 20:02:05 +0200 Subject: [PATCH 01/12] Make compute stress/mop compatible with 2D systems Issue an error if the stress should be computed in the Z direction in 2D systems The normalizing 'area' in 2D systems is the length of the simulation box in the other cartesian direction --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index ee8f5e554a..fb4725eeb4 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -148,7 +148,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( // 3D only - if (domain->dimension != 3) error->all(FLERR, "Compute stress/mop requires a 3d system"); + if (domain->dimension == 2 && dir == Z) + error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system"); // orthogonal simulation box if (domain->triclinic != 0) @@ -210,10 +211,14 @@ void ComputeStressMop::init() // Plane area - area = 1; - int i; - for (i = 0; i < 3; i++) { - if (i != dir) area = area * domain->prd[i]; + if (domain->dimension == 3) { + area = 1; + int i; + for (i = 0; i < 3; i++) { + if (i != dir) area = area * domain->prd[i]; + } + } else { + area = (dir == X) ? domain->prd[1] : domain->prd[0]; } // Timestep Value From e42aff54f906adadce6ac9e6fb6847e0e6ddb0b3 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 4 May 2024 20:07:08 +0200 Subject: [PATCH 02/12] Make compute stress/mop/profile compatible with 2D systems Issue an error if the stress is requested in the Z direction for 2D systems The normalizing 'area' is the length of the opposite cartesian direction --- .../compute_stress_mop_profile.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 676b0f5796..ca2d095fd9 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -133,8 +133,8 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a // 3D only - if (domain->dimension < 3) - error->all(FLERR, "Compute stress/mop/profile incompatible with simulation dimension"); + if (domain->dimension == 2 && dir == Z) + error->all(FLERR, "Compute stress/mop/profile incompatible with Z in 2d system"); // orthogonal simulation box @@ -198,11 +198,14 @@ void ComputeStressMopProfile::init() ftm2v = force->ftm2v; // plane area - - area = 1; - int i; - for (i = 0; i < 3; i++) { - if (i != dir) area = area * domain->prd[i]; + if (domain->dimension == 3) { + area = 1; + int i; + for (i = 0; i < 3; i++) { + if (i != dir) area = area * domain->prd[i]; + } + } else { + area = (dir == X) ? domain->prd[1] : domain->prd[0]; } // timestep Value From 46a441451ded609b102b9670a91c35312b3733aa Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 4 May 2024 20:11:08 +0200 Subject: [PATCH 03/12] Update compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 6630c7171f..521a017195 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -126,7 +126,7 @@ These styles are part of the EXTRA-COMPUTE package. They are only enabled if LAMMPS is built with that package. See the :doc:`Build package ` doc page on for more info. -The method is only implemented for 3d orthogonal simulation boxes whose +The method is only implemented for orthogonal simulation boxes whose size does not change in time, and axis-aligned planes. The method only works with two-body pair interactions, because it From 7b007d82a03bfdb5b1fa7d2b4fcfca8b803c70e2 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 9 May 2024 15:12:28 +0200 Subject: [PATCH 04/12] Make compute stress/mop compatible with triclinic boxes --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index fb4725eeb4..cd2069c7c3 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -89,12 +89,6 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( error->all(FLERR, "Plane for compute stress/mop is out of bounds"); } - if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) { - pos1 = pos + domain->prd[dir]; - } else { - pos1 = pos - domain->prd[dir]; - } - // parse values until one isn't recognized which = new int[3 * (narg - 5)]; @@ -151,10 +145,6 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( if (domain->dimension == 2 && dir == Z) error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system"); - // orthogonal simulation box - if (domain->triclinic != 0) - error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box"); - // Initialize some variables values_local = values_global = vector = nullptr; @@ -397,6 +387,9 @@ void ComputeStressMop::compute_pairs() jlist = firstneigh[i]; jnum = numneigh[i]; + xi[dir] -= pos; + domain->minimum_image(xi[0], xi[1], xi[2]); + for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; @@ -409,6 +402,10 @@ void ComputeStressMop::compute_pairs() xj[0] = atom->x[j][0]; xj[1] = atom->x[j][1]; xj[2] = atom->x[j][2]; + + xj[dir] -= pos; + domain->minimum_image(xj[0], xj[1], xj[2]); + delx = xi[0] - xj[0]; dely = xi[1] - xj[1]; delz = xi[2] - xj[2]; @@ -419,20 +416,19 @@ void ComputeStressMop::compute_pairs() if (newton_pair || j < nlocal) { //check if ij pair is across plane, add contribution to pressure - if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) { + if ((xi[dir] > 0.0) && (xj[dir] < 0.0)) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] += fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p; values_local[m + 2] += fpair * (xi[2] - xj[2]) / area * nktv2p; - } else if (((xi[dir] < pos) && (xj[dir] > pos)) || - ((xi[dir] < pos1) && (xj[dir] > pos1))) { + } else if ((xi[dir] < 0.0) && (xj[dir] > 0.0)) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] -= fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] -= fpair * (xi[1] - xj[1]) / area * nktv2p; values_local[m + 2] -= fpair * (xi[2] - xj[2]) / area * nktv2p; } } else { - if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) { + if ((xi[dir] > 0.0) && (xj[dir] < 0.0)) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] += fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p; From e2984c9724647c2df9bed1853e8444cc902f5bad Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 9 May 2024 15:13:26 +0200 Subject: [PATCH 05/12] Delete pos1 variable from compute_stress_mop.h --- src/EXTRA-COMPUTE/compute_stress_mop.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index 0a0ea8b55a..17c72826ab 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -51,7 +51,7 @@ class ComputeStressMop : public Compute { double *bond_local, *bond_global; double *angle_local, *angle_global; double *dihedral_local, *dihedral_global; - double pos, pos1, dt, nktv2p, ftm2v; + double pos, dt, nktv2p, ftm2v; double area; class NeighList *list; }; From 9aefa047cb0d1cc08f5976d861974421ee77ddfa Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 9 May 2024 15:17:30 +0200 Subject: [PATCH 06/12] Update compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 521a017195..6ee7e92720 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -126,8 +126,9 @@ These styles are part of the EXTRA-COMPUTE package. They are only enabled if LAMMPS is built with that package. See the :doc:`Build package ` doc page on for more info. -The method is only implemented for orthogonal simulation boxes whose -size does not change in time, and axis-aligned planes. +The method is implemented for simulation boxes whose +size does not change in time, and axis-aligned planes. Additionally, for +compute *stress/mop/profile*, the simulation box must be orthogonal. The method only works with two-body pair interactions, because it requires the class method ``Pair::single()`` to be implemented, which is From ac0513b5c4cb8a5a32e1afb20096246216219074 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 9 May 2024 15:25:57 +0200 Subject: [PATCH 07/12] whitespace in compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 6ee7e92720..e4127541a9 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -127,7 +127,7 @@ enabled if LAMMPS is built with that package. See the :doc:`Build package ` doc page on for more info. The method is implemented for simulation boxes whose -size does not change in time, and axis-aligned planes. Additionally, for +size does not change in time, and axis-aligned planes. Additionally, for compute *stress/mop/profile*, the simulation box must be orthogonal. The method only works with two-body pair interactions, because it From ed507dc676f9a7e7106946911d82ac55ca89667d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 13 May 2024 21:20:32 +0200 Subject: [PATCH 08/12] revert changes in compute_stress_mop.h --- src/EXTRA-COMPUTE/compute_stress_mop.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index 17c72826ab..0a0ea8b55a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -51,7 +51,7 @@ class ComputeStressMop : public Compute { double *bond_local, *bond_global; double *angle_local, *angle_global; double *dihedral_local, *dihedral_global; - double pos, dt, nktv2p, ftm2v; + double pos, pos1, dt, nktv2p, ftm2v; double area; class NeighList *list; }; From 86103fa89badfe304a07d7186cce5336dab7c7ad Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 13 May 2024 21:23:50 +0200 Subject: [PATCH 09/12] revert changes for triclinic boxes in compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index cd2069c7c3..32afa9975b 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -89,6 +89,12 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( error->all(FLERR, "Plane for compute stress/mop is out of bounds"); } + if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) { + pos1 = pos + domain->prd[dir]; + } else { + pos1 = pos - domain->prd[dir]; + } + // parse values until one isn't recognized which = new int[3 * (narg - 5)]; @@ -387,9 +393,6 @@ void ComputeStressMop::compute_pairs() jlist = firstneigh[i]; jnum = numneigh[i]; - xi[dir] -= pos; - domain->minimum_image(xi[0], xi[1], xi[2]); - for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; @@ -403,9 +406,6 @@ void ComputeStressMop::compute_pairs() xj[1] = atom->x[j][1]; xj[2] = atom->x[j][2]; - xj[dir] -= pos; - domain->minimum_image(xj[0], xj[1], xj[2]); - delx = xi[0] - xj[0]; dely = xi[1] - xj[1]; delz = xi[2] - xj[2]; @@ -416,19 +416,20 @@ void ComputeStressMop::compute_pairs() if (newton_pair || j < nlocal) { //check if ij pair is across plane, add contribution to pressure - if ((xi[dir] > 0.0) && (xj[dir] < 0.0)) { + if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] += fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p; values_local[m + 2] += fpair * (xi[2] - xj[2]) / area * nktv2p; - } else if ((xi[dir] < 0.0) && (xj[dir] > 0.0)) { + } else if (((xi[dir] < pos) && (xj[dir] > pos)) || + ((xi[dir] < pos1) && (xj[dir] > pos1))) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] -= fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] -= fpair * (xi[1] - xj[1]) / area * nktv2p; values_local[m + 2] -= fpair * (xi[2] - xj[2]) / area * nktv2p; } } else { - if ((xi[dir] > 0.0) && (xj[dir] < 0.0)) { + if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] += fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p; From 47c86cdf65fe5e48232e1d9eb1e6cf65784a2183 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 13 May 2024 21:26:01 +0200 Subject: [PATCH 10/12] Update compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index e4127541a9..b4779b8bf3 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -126,9 +126,8 @@ These styles are part of the EXTRA-COMPUTE package. They are only enabled if LAMMPS is built with that package. See the :doc:`Build package ` doc page on for more info. -The method is implemented for simulation boxes whose -size does not change in time, and axis-aligned planes. Additionally, for -compute *stress/mop/profile*, the simulation box must be orthogonal. +The method is implemented for orthogonal simulation boxes whose +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 From f05a551ffe370104132f4390be98f9cb78a6e548 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 13 May 2024 21:32:17 +0200 Subject: [PATCH 11/12] Update compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 32afa9975b..b8d21d2a4f 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -146,8 +146,11 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( // Error checks: - // 3D only + // orthogonal simulation box + if (domain->triclinic != 0) + error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box"); + // 2D and pressure calculation in the Z coordinate if (domain->dimension == 2 && dir == Z) error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system"); From 83a4ff06bd0e3af56dce6451e5977c4847c25120 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 14 May 2024 08:14:38 -0400 Subject: [PATCH 12/12] fix segfault --- src/REAXFF/fix_reaxff_species.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 2c0775b9d9..0183d2670b 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -346,13 +346,14 @@ int FixReaxFFSpecies::setmask() void FixReaxFFSpecies::setup(int /*vflag*/) { ntotal = static_cast(atom->natoms); - if (Name == nullptr) memory->create(Name, nutypes, "reaxff/species:Name"); if (!eleflag) { for (int i = 0; i < ntypes; i++) eletype[i] = reaxff->eletype[i+1]; GetUniqueElements(); } + memory->destroy(Name); + memory->create(Name, nutypes, "reaxff/species:Name"); post_integrate(); }