From 03d10e6bbca60d4e4f71909dd26fa4d00f2c4bcd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 18:56:07 -0400 Subject: [PATCH 01/24] correct documentation to add "inputs local" to compute reduce commands on local data --- doc/src/compute_angle_local.rst | 2 +- doc/src/compute_bond_local.rst | 2 +- doc/src/compute_dihedral_local.rst | 2 +- doc/src/fix_halt.rst | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/compute_angle_local.rst b/doc/src/compute_angle_local.rst index 126dd309af..d4491c6945 100644 --- a/doc/src/compute_angle_local.rst +++ b/doc/src/compute_angle_local.rst @@ -78,7 +78,7 @@ system and output the statistics in various ways: compute 2 all angle/local eng theta v_cos v_cossq set theta t dump 1 all local 100 tmp.dump c_1[*] c_2[*] - compute 3 all reduce ave c_2[*] + compute 3 all reduce ave c_2[*] inputs local thermo_style custom step temp press c_3[*] fix 10 all ave/histo 10 10 100 -1 1 20 c_2[3] mode vector file tmp.histo diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index 5036358c8c..e070d507b1 100644 --- a/doc/src/compute_bond_local.rst +++ b/doc/src/compute_bond_local.rst @@ -139,7 +139,7 @@ output the statistics in various ways: compute 2 all bond/local engpot dist v_dsq set dist d dump 1 all local 100 tmp.dump c_1[*] c_2[*] - compute 3 all reduce ave c_2[*] + compute 3 all reduce ave c_2[*] inputs local thermo_style custom step temp press c_3[*] fix 10 all ave/histo 10 10 100 0 6 20 c_2[3] mode vector file tmp.histo diff --git a/doc/src/compute_dihedral_local.rst b/doc/src/compute_dihedral_local.rst index 291b870373..d809cd39ce 100644 --- a/doc/src/compute_dihedral_local.rst +++ b/doc/src/compute_dihedral_local.rst @@ -76,7 +76,7 @@ angle in the system and output the statistics in various ways: compute 2 all dihedral/local phi v_cos v_cossq set phi p dump 1 all local 100 tmp.dump c_1[*] c_2[*] - compute 3 all reduce ave c_2[*] + compute 3 all reduce ave c_2[*] inputs local thermo_style custom step temp press c_3[*] fix 10 all ave/histo 10 10 100 -1 1 20 c_2[2] mode vector file tmp.histo diff --git a/doc/src/fix_halt.rst b/doc/src/fix_halt.rst index 25331804aa..0bcf2fb5ea 100644 --- a/doc/src/fix_halt.rst +++ b/doc/src/fix_halt.rst @@ -101,7 +101,7 @@ hstyle = bondmax option. .. code-block:: LAMMPS compute bdist all bond/local dist - compute bmax all reduce max c_bdist + compute bmax all reduce max c_bdist inputs local variable bondmax equal c_bmax Thus these two versions of a fix halt command will do the same thing: From b3b9d032fa1a83e3806b19f0f528c0be53f08ca5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 18:56:41 -0400 Subject: [PATCH 02/24] cleanup, clarification and re-wrap of doc file sections --- doc/src/Howto_chunk.rst | 67 ++++++++++++----------- doc/src/compute_reduce.rst | 10 ++-- doc/src/compute_reduce_chunk.rst | 91 +++++++++++++++++--------------- doc/src/group.rst | 4 +- 4 files changed, 93 insertions(+), 79 deletions(-) diff --git a/doc/src/Howto_chunk.rst b/doc/src/Howto_chunk.rst index f8655b745d..ea000eb22f 100644 --- a/doc/src/Howto_chunk.rst +++ b/doc/src/Howto_chunk.rst @@ -58,28 +58,30 @@ chunk ID for an individual atom can also be static (e.g. a molecule ID), or dynamic (e.g. what spatial bin an atom is in as it moves). Note that this compute allows the per-atom output of other -:doc:`computes `, :doc:`fixes `, and -:doc:`variables ` to be used to define chunk IDs for each -atom. This means you can write your own compute or fix to output a -per-atom quantity to use as chunk ID. See the :doc:`Modify ` -doc pages for info on how to do this. You can also define a :doc:`per-atom variable ` in the input script that uses a formula to -generate a chunk ID for each atom. +:doc:`computes `, :doc:`fixes `, and :doc:`variables +` to be used to define chunk IDs for each atom. This means +you can write your own compute or fix to output a per-atom quantity to +use as chunk ID. See the :doc:`Modify ` doc pages for info on +how to do this. You can also define a :doc:`per-atom variable +` in the input script that uses a formula to generate a chunk +ID for each atom. Fix ave/chunk command: ---------------------- -This fix takes the ID of a :doc:`compute chunk/atom ` command as input. For each chunk, -it then sums one or more specified per-atom values over the atoms in -each chunk. The per-atom values can be any atom property, such as -velocity, force, charge, potential energy, kinetic energy, stress, -etc. Additional keywords are defined for per-chunk properties like -density and temperature. More generally any per-atom value generated -by other :doc:`computes `, :doc:`fixes `, and :doc:`per-atom variables `, can be summed over atoms in each chunk. +This fix takes the ID of a :doc:`compute chunk/atom +` command as input. For each chunk, it then sums +one or more specified per-atom values over the atoms in each chunk. The +per-atom values can be any atom property, such as velocity, force, +charge, potential energy, kinetic energy, stress, etc. Additional +keywords are defined for per-chunk properties like density and +temperature. More generally any per-atom value generated by other +:doc:`computes `, :doc:`fixes `, and :doc:`per-atom +variables `, can be summed over atoms in each chunk. Similar to other averaging fixes, this fix allows the summed per-chunk values to be time-averaged in various ways, and output to a file. The -fix produces a global array as output with one row of values per -chunk. +fix produces a global array as output with one row of values per chunk. Compute \*/chunk commands: -------------------------- @@ -97,17 +99,20 @@ category: * :doc:`compute torque/chunk ` * :doc:`compute vcm/chunk ` -They each take the ID of a :doc:`compute chunk/atom ` command as input. As their names -indicate, they calculate the center-of-mass, radius of gyration, -moments of inertia, mean-squared displacement, temperature, torque, -and velocity of center-of-mass for each chunk of atoms. The :doc:`compute property/chunk ` command can tally the -count of atoms in each chunk and extract other per-chunk properties. +They each take the ID of a :doc:`compute chunk/atom +` command as input. As their names indicate, they +calculate the center-of-mass, radius of gyration, moments of inertia, +mean-squared displacement, temperature, torque, and velocity of +center-of-mass for each chunk of atoms. The :doc:`compute +property/chunk ` command can tally the count of +atoms in each chunk and extract other per-chunk properties. -The reason these various calculations are not part of the :doc:`fix ave/chunk command `, is that each requires a more +The reason these various calculations are not part of the :doc:`fix +ave/chunk command `, is that each requires a more complicated operation than simply summing and averaging over per-atom -values in each chunk. For example, many of them require calculation -of a center of mass, which requires summing mass\*position over the -atoms and then dividing by summed mass. +values in each chunk. For example, many of them require calculation of +a center of mass, which requires summing mass\*position over the atoms +and then dividing by summed mass. All of these computes produce a global vector or global array as output, with one or more values per chunk. The output can be used in @@ -118,9 +123,10 @@ various ways: * As input to the :doc:`fix ave/histo ` command to histogram values across chunks. E.g. a histogram of cluster sizes or molecule diffusion rates. -* As input to special functions of :doc:`equal-style variables `, like sum() and max() and ave(). E.g. to - find the largest cluster or fastest diffusing molecule or average - radius-of-gyration of a set of molecules (chunks). +* As input to special functions of :doc:`equal-style variables + `, like sum() and max() and ave(). E.g. to find the largest + cluster or fastest diffusing molecule or average radius-of-gyration of + a set of molecules (chunks). Other chunk commands: --------------------- @@ -138,9 +144,10 @@ spatially average per-chunk values calculated by a per-chunk compute. The :doc:`compute reduce/chunk ` command reduces a peratom value across the atoms in each chunk to produce a value per -chunk. When used with the :doc:`compute chunk/spread/atom ` command it can -create peratom values that induce a new set of chunks with a second -:doc:`compute chunk/atom ` command. +chunk. When used with the :doc:`compute chunk/spread/atom +` command it can create peratom values that +induce a new set of chunks with a second :doc:`compute chunk/atom +` command. Example calculations with chunks -------------------------------- diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index 9eafe7cd5a..f25a994764 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -206,11 +206,13 @@ IDs and the bond stretch will be printed with thermodynamic output. The *inputs* keyword allows selection of whether all the inputs are per-atom or local quantities. As noted above, all the inputs must be -the same kind (per-atom or local). Per-atom is the default setting. -If a compute or fix is specified as an input, it must produce per-atom -or local data to match this setting. If it produces both, e.g. for +the same kind (per-atom or local). Per-atom is the default setting. If +a compute or fix is specified as an input, it must produce per-atom or +local data to match this setting. If it produces both, like for example the :doc:`compute voronoi/atom ` command, then -this keyword selects between them. +this keyword selects between them. If a compute *only* produces local +data, like for example the :doc:`compute bond/local +command`, the setting "inputs local" is *required*. ---------- diff --git a/doc/src/compute_reduce_chunk.rst b/doc/src/compute_reduce_chunk.rst index 8ec19ade66..091604c0b2 100644 --- a/doc/src/compute_reduce_chunk.rst +++ b/doc/src/compute_reduce_chunk.rst @@ -37,55 +37,57 @@ Description Define a calculation that reduces one or more per-atom vectors into per-chunk values. This can be useful for diagnostic output. Or when -used in conjunction with the :doc:`compute chunk/spread/atom ` command it can be -used to create per-atom values that induce a new set of chunks with a -second :doc:`compute chunk/atom ` command. An -example is given below. +used in conjunction with the :doc:`compute chunk/spread/atom +` command it can be used to create per-atom +values that induce a new set of chunks with a second :doc:`compute +chunk/atom ` command. An example is given below. -In LAMMPS, chunks are collections of atoms defined by a :doc:`compute chunk/atom ` command, which assigns each atom -to a single chunk (or no chunk). The ID for this command is specified -as chunkID. For example, a single chunk could be the atoms in a -molecule or atoms in a spatial bin. See the :doc:`compute chunk/atom ` and :doc:`Howto chunk ` -doc pages for details of how chunks can be defined and examples of how -they can be used to measure properties of a system. +In LAMMPS, chunks are collections of atoms defined by a :doc:`compute +chunk/atom ` command, which assigns each atom to a +single chunk (or no chunk). The ID for this command is specified as +chunkID. For example, a single chunk could be the atoms in a molecule +or atoms in a spatial bin. See the :doc:`compute chunk/atom +` and :doc:`Howto chunk ` doc pages for +details of how chunks can be defined and examples of how they can be +used to measure properties of a system. For each atom, this compute accesses its chunk ID from the specified -*chunkID* compute. The per-atom value from an input contributes -to a per-chunk value corresponding the the chunk ID. +*chunkID* compute. The per-atom value from an input contributes to a +per-chunk value corresponding the the chunk ID. The reduction operation is specified by the *mode* setting and is performed over all the per-atom values from the atoms in each chunk. The *sum* option adds the pre-atom values to a per-chunk total. The -*min* or *max* options find the minimum or maximum value of the -per-atom values for each chunk. +*min* or *max* options find the minimum or maximum value of the per-atom +values for each chunk. -Note that only atoms in the specified group contribute to the -reduction operation. If the *chunkID* compute returns a 0 for the -chunk ID of an atom (i.e., the atom is not in a chunk defined by the -:doc:`compute chunk/atom ` command), that atom will -also not contribute to the reduction operation. An input that is a -compute or fix may define its own group which affects the quantities -it returns. For example, a compute with return a zero value for atoms -that are not in the group specified for that compute. +Note that only atoms in the specified group contribute to the reduction +operation. If the *chunkID* compute returns a 0 for the chunk ID of an +atom (i.e., the atom is not in a chunk defined by the :doc:`compute +chunk/atom ` command), that atom will also not +contribute to the reduction operation. An input that is a compute or +fix may define its own group which affects the quantities it returns. +For example, a compute with return a zero value for atoms that are not +in the group specified for that compute. Each listed input is operated on independently. Each input can be the -result of a :doc:`compute ` or :doc:`fix ` or the evaluation -of an atom-style :doc:`variable `. +result of a :doc:`compute ` or :doc:`fix ` or the +evaluation of an atom-style :doc:`variable `. -Note that for values from a compute or fix, the bracketed index I can -be specified using a wildcard asterisk with the index to effectively +Note that for values from a compute or fix, the bracketed index I can be +specified using a wildcard asterisk with the index to effectively specify multiple values. This takes the form "\*" or "\*n" or "m\*" or -"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or the -number of columns in the array (for *mode* = vector), then an asterisk -with no numeric values means all indices from 1 to :math:`N`. A leading -asterisk means all indices from 1 to n (inclusive). A trailing -asterisk means all indices from n to :math:`N` (inclusive). A middle asterisk -means all indices from m to n (inclusive). +"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or +the number of columns in the array (for *mode* = vector), then an +asterisk with no numeric values means all indices from 1 to :math:`N`. +A leading asterisk means all indices from 1 to n (inclusive). A +trailing asterisk means all indices from n to :math:`N` (inclusive). A +middle asterisk means all indices from m to n (inclusive). Using a wildcard is the same as if the individual columns of the array -had been listed one by one. For example, the following two compute reduce/chunk -commands are equivalent, since the -:doc:`compute property/chunk ` command creates a per-atom +had been listed one by one. For example, the following two compute +reduce/chunk commands are equivalent, since the :doc:`compute +property/chunk ` command creates a per-atom array with 3 columns: .. code-block:: LAMMPS @@ -164,13 +166,14 @@ Output info """"""""""" This compute calculates a global vector if a single input value is -specified, otherwise a global array is output. The number of columns -in the array is the number of inputs provided. The length of the -vector or the number of vector elements or array rows = the number of -chunks *Nchunk* as calculated by the specified :doc:`compute chunk/atom ` command. The vector or array can -be accessed by any command that uses global values from a compute as -input. See the :doc:`Howto output ` page for an -overview of LAMMPS output options. +specified, otherwise a global array is output. The number of columns in +the array is the number of inputs provided. The length of the vector or +the number of vector elements or array rows = the number of chunks +*Nchunk* as calculated by the specified :doc:`compute chunk/atom +` command. The vector or array can be accessed by +any command that uses global values from a compute as input. See the +:doc:`Howto output ` page for an overview of LAMMPS output +options. The per-atom values for the vector or each column of the array will be in whatever :doc:`units ` the corresponding input value is in. @@ -183,7 +186,9 @@ Restrictions Related commands """""""""""""""" -:doc:`compute chunk/atom `, :doc:`compute reduce `, :doc:`compute chunk/spread/atom ` +:doc:`compute chunk/atom `, +:doc:`compute reduce `, +:doc:`compute chunk/spread/atom ` Default """"""" diff --git a/doc/src/group.rst b/doc/src/group.rst index 15ab0c9dc8..a7a29467ff 100644 --- a/doc/src/group.rst +++ b/doc/src/group.rst @@ -162,7 +162,7 @@ potential energy is above the threshold value :math:`-3.0`. compute 1 all pe/atom compute 2 all reduce sum c_1 thermo_style custom step temp pe c_2 - run 0 + run 0 post no variable eatom atom "c_1 > -3.0" group hienergy variable eatom @@ -173,7 +173,7 @@ Note that these lines compute 2 all reduce sum c_1 thermo_style custom step temp pe c_2 - run 0 + run 0 post no are necessary to ensure that the "eatom" variable is current when the group command invokes it. Because the eatom variable computes the From eafef460e289f004230f3f1300507619ea1e206c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 19:24:35 -0400 Subject: [PATCH 03/24] silence compiler warning --- src/DIELECTRIC/pppm_disp_dielectric.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp index 346df8bf93..aba8ccf486 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.cpp +++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp @@ -715,7 +715,7 @@ void PPPMDispDielectric::fieldforce_c_ad() { int i,l,m,n,nx,ny,nz,mx,my,mz; FFT_SCALAR dx,dy,dz; - FFT_SCALAR ekx,eky,ekz,u; + FFT_SCALAR ekx,eky,ekz; double s1,s2,s3; double sf = 0.0; From c15c4d23bb50bdfc19cc2a2a15e2a739a1de7c58 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 19:30:30 -0400 Subject: [PATCH 04/24] fix spelling/syntax issues --- doc/src/compute_reduce.rst | 4 ++-- doc/utils/sphinx-config/false_positives.txt | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index f25a994764..e5c99a478f 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -211,8 +211,8 @@ a compute or fix is specified as an input, it must produce per-atom or local data to match this setting. If it produces both, like for example the :doc:`compute voronoi/atom ` command, then this keyword selects between them. If a compute *only* produces local -data, like for example the :doc:`compute bond/local -command`, the setting "inputs local" is *required*. +data, like for example the :doc:`compute bond/local command +`, the setting "inputs local" is *required*. ---------- diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 2754e4bb1c..0708a14882 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1286,6 +1286,7 @@ gcc gcmc gdot GeC +Geesthacht Geier gencode Geocomputing @@ -4219,6 +4220,7 @@ zeeman Zeeman Zemer zenodo +Zentrum Zepeda zflag Zhang From 8156e5661749dfb3b6d4e82b92d4f909e2571365 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 19:31:40 -0400 Subject: [PATCH 05/24] be more paranoid about Sphinx updates breaking extensions we use --- doc/utils/requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/utils/requirements.txt b/doc/utils/requirements.txt index 4b9e1a4e48..acf575fe58 100644 --- a/doc/utils/requirements.txt +++ b/doc/utils/requirements.txt @@ -1,4 +1,4 @@ -Sphinx >= 5.3.0, <9.0 +Sphinx >= 5.3.0, <8.2.0 sphinxcontrib-spelling sphinxcontrib-jquery sphinx-design From f7291713f534c981c513a662c61e9c25a7e2e295 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 19:34:56 -0400 Subject: [PATCH 06/24] error out when fix nve is used with additional arguments --- src/fix_nve.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index 30c1f24867..58a5677b92 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -15,6 +15,7 @@ #include "fix_nve.h" #include "atom.h" +#include "error.h" #include "force.h" #include "respa.h" #include "update.h" @@ -28,7 +29,11 @@ FixNVE::FixNVE(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (!utils::strmatch(style,"^nve/sphere") && narg < 3) - utils::missing_cmd_args(FLERR, "fix nve", error); + utils::missing_cmd_args(FLERR, fmt::format("fix {}", style), error); + + auto plain_style = utils::strip_style_suffix(style, lmp); + if (utils::strmatch(plain_style, "^nve$") && narg > 3) + error->all(FLERR, "Unsupported additional arguments for fix {}", style); dynamic_group_allow = 1; time_integrate = 1; From 1c9daad65749bda99af7974db25e310965beb644 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 21:57:17 -0400 Subject: [PATCH 07/24] change tutorial download URL to tutorial website --- tools/lammps-gui/lammpsgui.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/lammps-gui/lammpsgui.cpp b/tools/lammps-gui/lammpsgui.cpp index fe6b8c5391..c8a715cf6b 100644 --- a/tools/lammps-gui/lammpsgui.cpp +++ b/tools/lammps-gui/lammpsgui.cpp @@ -2016,8 +2016,8 @@ bool LammpsGui::eventFilter(QObject *watched, QEvent *event) } // LAMMPS geturl command with current location of the input and solution files on the web -static const QString geturl = "geturl https://raw.githubusercontent.com/akohlmey/" - "lammps-tutorials-inputs/main/tutorial%1/%2 output %2 verify no"; +static const QString geturl = "geturl https://raw.githubusercontent.com/lammpstutorials/" + "lammpstutorials-article/refs/heads/main/files/tutorial%1/%2 output %2 verify no"; void LammpsGui::setup_tutorial(int tutno, const QString &dir, bool purgedir, bool getsolution) { From 1e63f031f0e7348cad2642a7e6c72eacec3d32d5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 13 Oct 2024 22:52:40 -0400 Subject: [PATCH 08/24] remove empirical filter to remove outliers from corrupted data --- tools/lammps-gui/chartviewer.cpp | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/tools/lammps-gui/chartviewer.cpp b/tools/lammps-gui/chartviewer.cpp index eb5444de8e..8543da166a 100644 --- a/tools/lammps-gui/chartviewer.cpp +++ b/tools/lammps-gui/chartviewer.cpp @@ -423,27 +423,6 @@ void ChartViewer::add_data(int step, double data) { if (last_step < step) { last_step = step; - - // do not add data that deviates by more than 4 sigma from the average - // over the last 5 to 20 data items. this is a hack to work around - // getting corrupted data from lammps_get_last_thermo() - const auto &points = series->points(); - const auto count = points.count(); - if (count > 4) { - double ysum = 0.0; - double ysumsq = 0.0; - int first = count - 20; - if (first < 0) first = 0; - for (int i = first; i < count; ++i) { - double val = points[i].y(); - ysum += val; - ysumsq += val * val; - } - const double num = count - first; - const double avg = ysum / num; - const double avgsq = ysumsq / num; - if (fabs(data - avg) > (5.0 * sqrt(avgsq - (avg * avg)))) return; - } series->append(step, data); QSettings settings; From db3416c4b30284a221937dd32fba69b6c6f6821c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 14 Oct 2024 08:19:53 -0400 Subject: [PATCH 09/24] Apply corrections from code review by @simongravelle Co-authored-by: Simon Gravelle --- doc/src/compute_reduce_chunk.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/compute_reduce_chunk.rst b/doc/src/compute_reduce_chunk.rst index 091604c0b2..eeadd50a97 100644 --- a/doc/src/compute_reduce_chunk.rst +++ b/doc/src/compute_reduce_chunk.rst @@ -53,11 +53,11 @@ used to measure properties of a system. For each atom, this compute accesses its chunk ID from the specified *chunkID* compute. The per-atom value from an input contributes to a -per-chunk value corresponding the the chunk ID. +per-chunk value corresponding the chunk ID. The reduction operation is specified by the *mode* setting and is performed over all the per-atom values from the atoms in each chunk. -The *sum* option adds the pre-atom values to a per-chunk total. The +The *sum* option adds the per-atom values to a per-chunk total. The *min* or *max* options find the minimum or maximum value of the per-atom values for each chunk. @@ -67,7 +67,7 @@ atom (i.e., the atom is not in a chunk defined by the :doc:`compute chunk/atom ` command), that atom will also not contribute to the reduction operation. An input that is a compute or fix may define its own group which affects the quantities it returns. -For example, a compute with return a zero value for atoms that are not +For example, a compute will return a zero value for atoms that are not in the group specified for that compute. Each listed input is operated on independently. Each input can be the From 38500c647b42bb0aa76f9bd91f02074636f05d6d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 14 Oct 2024 09:41:45 -0400 Subject: [PATCH 10/24] move cursor to end of log buffer before inserting new text --- tools/lammps-gui/lammpsgui.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tools/lammps-gui/lammpsgui.cpp b/tools/lammps-gui/lammpsgui.cpp index c8a715cf6b..fb38b6fc92 100644 --- a/tools/lammps-gui/lammpsgui.cpp +++ b/tools/lammps-gui/lammpsgui.cpp @@ -988,6 +988,7 @@ void LammpsGui::logupdate() if (logwindow) { const auto text = capturer->GetChunk(); if (text.size() > 0) { + logwindow->moveCursor(QTextCursor::End); logwindow->insertPlainText(text.c_str()); logwindow->moveCursor(QTextCursor::End); logwindow->textCursor().deleteChar(); From 2ce934fdb9df2fed81512343da775f7a13b5b388 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 16 Oct 2024 10:05:26 -0400 Subject: [PATCH 11/24] added section on OPLS-AA taken from https://matsci.org/t/issues-with-opls-force-field-for-ethylene-glycol-simulation/58296/11 --- doc/src/Howto_bioFF.rst | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 72dcec2d31..94209d3515 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -1,5 +1,5 @@ -CHARMM, AMBER, COMPASS, and DREIDING force fields -================================================= +CHARMM, AMBER, COMPASS, DREIDING, and OPLS-AA force fields +========================================================== A compact summary of the concepts, definitions, and properties of force fields with explicit bonded interactions (like the ones discussed @@ -236,6 +236,36 @@ documentation for the formula it computes. * :doc:`special_bonds ` dreiding +OPLS-AA +------- + +OPLS-AA is a general force field for atomistic simulation of proteins in +their native environment, developed at the `Jorgensen group +`_ at Purdue University and +later at Yale University. + +This force field is based on atom types mapping to specific functional +groups in organic and biological molecules. The atom type includes an +average atomic charge reflecting the oxidation state of the element in a +specific chemical bond :ref:`(Jorgensen) `. + +The interaction styles listed below compute force field formulas that +are consistent with the OPLS-AA force field. See each command's +documentation for the formula it computes. + +* :doc:`bond_style ` harmonic +* :doc:`angle_style ` harmonic +* :doc:`dihedral_style ` opls +* :doc:`improper_style ` harmonic +* :doc:`pair_style ` buck +* :doc:`pair_style ` buck/coul/cut +* :doc:`pair_style ` buck/coul/long +* :doc:`pair_style ` lj/cut +* :doc:`pair_style ` lj/cut/coul/cut +* :doc:`pair_style ` lj/cut/coul/long +* :doc:`pair_modify ` geometric +* :doc:`special_bonds ` lj/coul 0.0 0.0 0.5 + ---------- .. _Typelabel2: @@ -266,3 +296,6 @@ documentation for the formula it computes. **(Mayo)** Mayo, Olfason, Goddard III (1990). J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010 +.. _howto-Jorgensen: + +**(Jorgensen)** Jorgensen, Tirado-Rives (1988). J Am Chem Soc, 110, 1657-1666. https://doi.org/10.1021/ja00214a001 From 86e18b8a249a6c2a77b09c4668ddf6c24d56952d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 16 Oct 2024 10:23:16 -0400 Subject: [PATCH 12/24] fix copy-n-modify issue --- doc/src/Howto_bioFF.rst | 4 ---- 1 file changed, 4 deletions(-) diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 94209d3515..23a235f454 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -257,10 +257,6 @@ documentation for the formula it computes. * :doc:`angle_style ` harmonic * :doc:`dihedral_style ` opls * :doc:`improper_style ` harmonic -* :doc:`pair_style ` buck -* :doc:`pair_style ` buck/coul/cut -* :doc:`pair_style ` buck/coul/long -* :doc:`pair_style ` lj/cut * :doc:`pair_style ` lj/cut/coul/cut * :doc:`pair_style ` lj/cut/coul/long * :doc:`pair_modify ` geometric From 2074b0c400c98c9b8d31ed712712068f43816a0c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Oct 2024 21:49:55 -0400 Subject: [PATCH 13/24] add check to allow only valid swaps between hybrid sub-styles and document it --- doc/src/fix_atom_swap.rst | 13 ++++++++----- src/MC/fix_atom_swap.cpp | 23 +++++++++++++++++++++++ src/pair_hybrid.h | 1 + 3 files changed, 32 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_atom_swap.rst b/doc/src/fix_atom_swap.rst index 48f5b33a74..1a833f69b1 100644 --- a/doc/src/fix_atom_swap.rst +++ b/doc/src/fix_atom_swap.rst @@ -119,15 +119,14 @@ groups of atoms that have different charges, these charges will not be changed when the atom types change. Since this fix computes total potential energies before and after -proposed swaps, so even complicated potential energy calculations are -OK, including the following: +proposed swaps, even complicated potential energy calculations are +acceptable, including the following: * long-range electrostatics (:math:`k`-space) * many body pair styles -* hybrid pair styles -* eam pair styles +* hybrid pair styles (with restrictions) +* EAM pair styles * triclinic systems -* need to include potential energy contributions from other fixes Some fixes have an associated potential energy. Examples of such fixes include: :doc:`efield `, :doc:`gravity `, @@ -181,6 +180,10 @@ This fix is part of the MC package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` doc page for more info. +When this fix is used with a :doc:`hybrid pair style ` +system, only swaps between atom types of the same substyle (or +combination of substyles) are permitted. + This fix cannot be used with systems that do not have per-type masses (e.g. atom style sphere) since the implemented algorithm pre-computes velocity rescaling factors from per-type masses and ignores any per-atom diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp index 1348f1248e..209c5cf8c4 100644 --- a/src/MC/fix_atom_swap.cpp +++ b/src/MC/fix_atom_swap.cpp @@ -35,6 +35,7 @@ #include "modify.h" #include "neighbor.h" #include "pair.h" +#include "pair_hybrid.h" #include "random_park.h" #include "region.h" #include "update.h" @@ -223,6 +224,28 @@ void FixAtomSwap::init() error->all(FLERR, "Mu not allowed when not using semi-grand in fix atom/swap command"); } + // check if constraints for hybrid pair styles are fulfilled + + if (utils::strmatch(force->pair_style, "^hybrid")) { + auto *hybrid = dynamic_cast(force->pair); + if (hybrid) { + for (int i = 0; i < nswaptypes - 1; ++i) { + int type1 = type_list[i]; + for (int j = i + 1; j < nswaptypes; ++j) { + int type2 = type_list[j]; + if (hybrid->nmap[type1][type1] != hybrid->nmap[type2][type2]) + error->all(FLERR, "Pair {} substyles for atom types {} and {} are not compatible", + force->pair_style, type1, type2); + for (int k = 0; k < hybrid->nmap[type1][type1]; ++k) { + if (hybrid->map[type1][type1][k] != hybrid->map[type2][type2][k]) + error->all(FLERR, "Pair {} substyles for atom types {} and {} are not compatible", + force->pair_style, type1, type2); + } + } + } + } + } + // set index and check validity of region if (idregion) { diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index b59cc82cb4..995de37e59 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -28,6 +28,7 @@ namespace LAMMPS_NS { class PairHybrid : public Pair { friend class AtomVecDielectric; friend class ComputeSpin; + friend class FixAtomSwap; friend class FixGPU; friend class FixIntel; friend class FixNVESpin; From f128abe00239d5be478847a4b732e1fecf206f78 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Oct 2024 21:51:43 -0400 Subject: [PATCH 14/24] improve OPLS description --- doc/src/Howto_bioFF.rst | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 23a235f454..31c509de7a 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -236,26 +236,33 @@ documentation for the formula it computes. * :doc:`special_bonds ` dreiding -OPLS-AA -------- +OPLS +---- -OPLS-AA is a general force field for atomistic simulation of proteins in -their native environment, developed at the `Jorgensen group +OPLS (Optimized Potentials for Liquid Simulations) is a general force +field for atomistic simulation of organic molecules in solvent. It was +developed at the `Jorgensen group `_ at Purdue University and -later at Yale University. +later at Yale University. Multiple versions of the OPLS parameters +exist for united atom representations (OPLS-UA) and for all-atom +representations (OPLS-AA). This force field is based on atom types mapping to specific functional groups in organic and biological molecules. The atom type includes an average atomic charge reflecting the oxidation state of the element in a -specific chemical bond :ref:`(Jorgensen) `. +specific chemical bond :ref:`(Jorgensen) ` and computed +based on increments determined by the atom type of the atoms bond to it. The interaction styles listed below compute force field formulas that -are consistent with the OPLS-AA force field. See each command's -documentation for the formula it computes. +are fully or in part consistent with the OPLS style force fields. See +each command's documentation for the formula it computes. Some are only +compatible with a subset of OPLS interactions. * :doc:`bond_style ` harmonic * :doc:`angle_style ` harmonic * :doc:`dihedral_style ` opls +* :doc:`improper_style ` cvff +* :doc:`improper_style ` fourier * :doc:`improper_style ` harmonic * :doc:`pair_style ` lj/cut/coul/cut * :doc:`pair_style ` lj/cut/coul/long From 34ab2f862a38790f9e5571bae93fb1c48d980232 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Oct 2024 22:09:25 -0400 Subject: [PATCH 15/24] small tweaks --- doc/src/Howto_bioFF.rst | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 31c509de7a..d35913eb0d 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -1,5 +1,5 @@ -CHARMM, AMBER, COMPASS, DREIDING, and OPLS-AA force fields -========================================================== +CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields +======================================================= A compact summary of the concepts, definitions, and properties of force fields with explicit bonded interactions (like the ones discussed @@ -241,17 +241,18 @@ OPLS OPLS (Optimized Potentials for Liquid Simulations) is a general force field for atomistic simulation of organic molecules in solvent. It was -developed at the `Jorgensen group +developed by the `Jorgensen group `_ at Purdue University and later at Yale University. Multiple versions of the OPLS parameters exist for united atom representations (OPLS-UA) and for all-atom representations (OPLS-AA). -This force field is based on atom types mapping to specific functional -groups in organic and biological molecules. The atom type includes an -average atomic charge reflecting the oxidation state of the element in a -specific chemical bond :ref:`(Jorgensen) ` and computed -based on increments determined by the atom type of the atoms bond to it. +This force field is based on atom types mapped to specific functional +groups in organic and biological molecules. Each atom includes an +static, partial atomic charge reflecting the oxidation state of the +element derived from its bonded neighbors :ref:`(Jorgensen) ` +and computed based on increments determined by the atom type of the +atoms bond to it. The interaction styles listed below compute force field formulas that are fully or in part consistent with the OPLS style force fields. See From 760d871b7ad96f78dcb9ed38235dc95157440a10 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 19 Oct 2024 11:02:53 -0400 Subject: [PATCH 16/24] enabled and apply clang-format --- src/compute_angle_local.cpp | 156 ++++++++------- src/compute_bond_local.cpp | 344 ++++++++++++++++++--------------- src/compute_dihedral_local.cpp | 176 +++++++++-------- src/compute_improper_local.cpp | 89 +++++---- 4 files changed, 405 insertions(+), 360 deletions(-) diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index 3e8b15fd64..569c30772e 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,38 +12,40 @@ ------------------------------------------------------------------------- */ #include "compute_angle_local.h" -#include -#include + +#include "angle.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "update.h" #include "domain.h" +#include "error.h" #include "force.h" -#include "angle.h" #include "input.h" -#include "variable.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "molecule.h" +#include "update.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; static constexpr int DELTA = 10000; -enum{THETA,ENG,VARIABLE}; +enum { THETA, ENG, VARIABLE }; /* ---------------------------------------------------------------------- */ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), + vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute angle/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute angle/local command"); if (atom->avec->angles_allow == 0) - error->all(FLERR,"Compute angle/local used when angles are not allowed"); + error->all(FLERR, "Compute angle/local used when angles are not allowed"); local_flag = 1; @@ -52,7 +53,7 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = narg - 3; bstyle = new int[nvalues]; - vstr = new char*[nvalues]; + vstr = new char *[nvalues]; vvar = new int[nvalues]; nvalues = 0; @@ -61,16 +62,17 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"theta") == 0) { + if (strcmp(arg[iarg], "theta") == 0) { bstyle[nvalues++] = THETA; tflag = 1; - } else if (strcmp(arg[iarg],"eng") == 0) { + } else if (strcmp(arg[iarg], "eng") == 0) { bstyle[nvalues++] = ENG; - } else if (strncmp(arg[iarg],"v_",2) == 0) { + } else if (strncmp(arg[iarg], "v_", 2) == 0) { bstyle[nvalues++] = VARIABLE; vstr[nvar] = utils::strdup(&arg[iarg][2]); nvar++; - } else break; + } else + break; } // optional args @@ -79,45 +81,46 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : tstr = nullptr; while (iarg < narg) { - if (strcmp(arg[iarg],"set") == 0) { + if (strcmp(arg[iarg], "set") == 0) { setflag = 1; - if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command"); - if (strcmp(arg[iarg+1],"theta") == 0) { - delete [] tstr; - tstr = utils::strdup(arg[iarg+2]); + if (iarg + 3 > narg) error->all(FLERR, "Illegal compute angle/local command"); + if (strcmp(arg[iarg + 1], "theta") == 0) { + delete[] tstr; + tstr = utils::strdup(arg[iarg + 2]); tflag = 1; - } else error->all(FLERR,"Illegal compute angle/local command"); + } else + error->all(FLERR, "Illegal compute angle/local command"); iarg += 3; - } else error->all(FLERR,"Illegal compute angle/local command"); + } else + error->all(FLERR, "Illegal compute angle/local command"); } // error check if (nvar) { - if (!setflag) - error->all(FLERR,"Compute angle/local variable requires a set variable"); + if (!setflag) error->all(FLERR, "Compute angle/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR,"Variable name for copute angle/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for copute angle/local does not exist"); if (!input->variable->equalstyle(vvar[i])) - error->all(FLERR,"Variable for compute angle/local is invalid style"); + error->all(FLERR, "Variable for compute angle/local is invalid style"); } if (tstr) { tvar = input->variable->find(tstr); - if (tvar < 0) - error->all(FLERR,"Variable name for compute angle/local does not exist"); + if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); if (!input->variable->internalstyle(tvar)) - error->all(FLERR,"Variable for compute angle/local is invalid style"); + error->all(FLERR, "Variable for compute angle/local is invalid style"); } } else if (setflag) - error->all(FLERR,"Compute angle/local set with no variable"); + error->all(FLERR, "Compute angle/local set with no variable"); // initialize output - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -128,12 +131,12 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : ComputeAngleLocal::~ComputeAngleLocal() { - delete [] bstyle; - for (int i = 0; i < nvar; i++) delete [] vstr[i]; - delete [] vstr; - delete [] vvar; + delete[] bstyle; + for (int i = 0; i < nvar; i++) delete[] vstr[i]; + delete[] vstr; + delete[] vvar; - delete [] tstr; + delete[] tstr; memory->destroy(vlocal); memory->destroy(alocal); @@ -144,19 +147,17 @@ ComputeAngleLocal::~ComputeAngleLocal() void ComputeAngleLocal::init() { if (force->angle == nullptr) - error->all(FLERR,"No angle style is defined for compute angle/local"); + error->all(FLERR, "No angle style is defined for compute angle/local"); if (nvar) { for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR,"Variable name for compute angle/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); } if (tstr) { tvar = input->variable->find(tstr); - if (tvar < 0) - error->all(FLERR,"Variable name for compute angle/local does not exist"); + if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); } } @@ -194,10 +195,10 @@ void ComputeAngleLocal::compute_local() int ComputeAngleLocal::compute_angles(int flag) { - int i,m,na,atom1,atom2,atom3,imol,iatom,atype,ivar; + int i, m, na, atom1, atom2, atom3, imol, iatom, atype, ivar; tagint tagprev; - double delx1,dely1,delz1,delx2,dely2,delz2; - double rsq1,rsq2,r1,r2,c,theta; + double delx1, dely1, delz1, delx2, dely2, delz2; + double rsq1, rsq2, r1, r2, c, theta; double *ptr; double **x = atom->x; @@ -224,7 +225,8 @@ int ComputeAngleLocal::compute_angles(int flag) for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) na = num_angle[atom2]; + if (molecular == Atom::MOLECULAR) + na = num_angle[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -242,8 +244,8 @@ int ComputeAngleLocal::compute_angles(int flag) 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); + 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; @@ -261,50 +263,54 @@ int ComputeAngleLocal::compute_angles(int flag) delx1 = x[atom1][0] - x[atom2][0]; dely1 = x[atom1][1] - x[atom2][1]; delz1 = x[atom1][2] - x[atom2][2]; - domain->minimum_image(delx1,dely1,delz1); + domain->minimum_image(delx1, dely1, delz1); - rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1; r1 = sqrt(rsq1); delx2 = x[atom3][0] - x[atom2][0]; dely2 = x[atom3][1] - x[atom2][1]; delz2 = x[atom3][2] - x[atom2][2]; - domain->minimum_image(delx2,dely2,delz2); + domain->minimum_image(delx2, dely2, delz2); - rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2; r2 = sqrt(rsq2); // c = cosine of angle // theta = angle in radians - c = delx1*delx2 + dely1*dely2 + delz1*delz2; - c /= r1*r2; + c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + c /= r1 * r2; if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; theta = acos(c); } - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; if (nvar) { ivar = 0; - if (tstr) input->variable->internal_set(tvar,theta); + if (tstr) input->variable->internal_set(tvar, theta); } for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { - case THETA: - ptr[n] = 180.0*theta/MY_PI; - break; - case ENG: - if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3); - else ptr[n] = 0.0; - break; - case VARIABLE: - ptr[n] = input->variable->compute_equal(vvar[ivar]); - ivar++; - break; + case THETA: + ptr[n] = 180.0 * theta / MY_PI; + break; + case ENG: + if (atype > 0) + ptr[n] = angle->single(atype, atom1, atom2, atom3); + else + ptr[n] = 0.0; + break; + case VARIABLE: + ptr[n] = input->variable->compute_equal(vvar[ivar]); + ivar++; + break; } } @@ -325,11 +331,11 @@ void ComputeAngleLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"angle/local:vector_local"); + memory->create(vlocal, nmax, "angle/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"angle/local:array_local"); + memory->create(alocal, nmax, nvalues, "angle/local:array_local"); array_local = alocal; } } @@ -340,6 +346,6 @@ void ComputeAngleLocal::reallocate(int n) double ComputeAngleLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index e3c3e26a85..bd59d76b9c 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -35,18 +34,35 @@ using namespace LAMMPS_NS; static constexpr int DELTA = 10000; -enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE,BN}; +enum { + DIST, + DX, + DY, + DZ, + VELVIB, + OMEGA, + ENGTRANS, + ENGVIB, + ENGROT, + ENGPOT, + FORCE, + FX, + FY, + FZ, + VARIABLE, + BN +}; /* ---------------------------------------------------------------------- */ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr), + vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute bond/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute bond/local command"); if (atom->avec->bonds_allow == 0) - error->all(FLERR,"Compute bond/local used when bonds are not allowed"); + error->all(FLERR, "Compute bond/local used when bonds are not allowed"); local_flag = 1; comm_forward = 3; @@ -56,7 +72,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = narg - 3; bstyle = new int[nvalues]; bindex = new int[nvalues]; - vstr = new char*[nvalues]; + vstr = new char *[nvalues]; vvar = new int[nvalues]; nvalues = 0; @@ -64,30 +80,45 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST; - else if (strcmp(arg[iarg],"dx") == 0) bstyle[nvalues++] = DX; - else if (strcmp(arg[iarg],"dy") == 0) bstyle[nvalues++] = DY; - else if (strcmp(arg[iarg],"dz") == 0) bstyle[nvalues++] = DZ; - else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT; - else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE; - else if (strcmp(arg[iarg],"fx") == 0) bstyle[nvalues++] = FX; - else if (strcmp(arg[iarg],"fy") == 0) bstyle[nvalues++] = FY; - else if (strcmp(arg[iarg],"fz") == 0) bstyle[nvalues++] = FZ; - else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB; - else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT; - else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS; - else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA; - else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB; - else if (strncmp(arg[iarg],"v_",2) == 0) { + if (strcmp(arg[iarg], "dist") == 0) + bstyle[nvalues++] = DIST; + else if (strcmp(arg[iarg], "dx") == 0) + bstyle[nvalues++] = DX; + else if (strcmp(arg[iarg], "dy") == 0) + bstyle[nvalues++] = DY; + else if (strcmp(arg[iarg], "dz") == 0) + bstyle[nvalues++] = DZ; + else if (strcmp(arg[iarg], "engpot") == 0) + bstyle[nvalues++] = ENGPOT; + else if (strcmp(arg[iarg], "force") == 0) + bstyle[nvalues++] = FORCE; + else if (strcmp(arg[iarg], "fx") == 0) + bstyle[nvalues++] = FX; + else if (strcmp(arg[iarg], "fy") == 0) + bstyle[nvalues++] = FY; + else if (strcmp(arg[iarg], "fz") == 0) + bstyle[nvalues++] = FZ; + else if (strcmp(arg[iarg], "engvib") == 0) + bstyle[nvalues++] = ENGVIB; + else if (strcmp(arg[iarg], "engrot") == 0) + bstyle[nvalues++] = ENGROT; + else if (strcmp(arg[iarg], "engtrans") == 0) + bstyle[nvalues++] = ENGTRANS; + else if (strcmp(arg[iarg], "omega") == 0) + bstyle[nvalues++] = OMEGA; + else if (strcmp(arg[iarg], "velvib") == 0) + bstyle[nvalues++] = VELVIB; + else if (strncmp(arg[iarg], "v_", 2) == 0) { bstyle[nvalues++] = VARIABLE; vstr[nvar] = utils::strdup(&arg[iarg][2]); nvar++; - } else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN + } else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN int n = std::stoi(&arg[iarg][1]); if (n <= 0) error->all(FLERR, "Invalid keyword {} in compute bond/local command", arg[iarg]); bstyle[nvalues] = BN; bindex[nvalues++] = n - 1; - } else break; + } else + break; } // optional args @@ -96,40 +127,39 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : dstr = nullptr; while (iarg < narg) { - if (strcmp(arg[iarg],"set") == 0) { + if (strcmp(arg[iarg], "set") == 0) { setflag = 1; - if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"compute bond/local set", error); - if (strcmp(arg[iarg+1],"dist") == 0) { - delete [] dstr; - dstr = utils::strdup(arg[iarg+2]); - } else error->all(FLERR,"Unknown compute bond/local set keyword: {}", arg[iarg+2]); + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "compute bond/local set", error); + if (strcmp(arg[iarg + 1], "dist") == 0) { + delete[] dstr; + dstr = utils::strdup(arg[iarg + 2]); + } else + error->all(FLERR, "Unknown compute bond/local set keyword: {}", arg[iarg + 2]); iarg += 3; - } else error->all(FLERR,"Unknown compute bond/local keyword: {}", arg[iarg]); + } else + error->all(FLERR, "Unknown compute bond/local keyword: {}", arg[iarg]); } // error check if (nvar) { - if (!setflag) - error->all(FLERR,"Compute bond/local variable requires a set variable"); + if (!setflag) error->all(FLERR, "Compute bond/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); if (vvar[i] < 0) - error->all(FLERR,"Variable name {} for copute bond/local does not exist", vstr[i]); + error->all(FLERR, "Variable name {} for copute bond/local does not exist", vstr[i]); if (!input->variable->equalstyle(vvar[i])) - error->all(FLERR,"Variable {} for compute bond/local is invalid style", vstr[i]); + error->all(FLERR, "Variable {} for compute bond/local is invalid style", vstr[i]); } if (dstr) { dvar = input->variable->find(dstr); - if (dvar < 0) - error->all(FLERR,"Variable name for compute bond/local does not exist"); + if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist"); if (!input->variable->internalstyle(dvar)) - error->all(FLERR,"Variable for compute bond/local is invalid style"); + error->all(FLERR, "Variable for compute bond/local is invalid style"); } } else if (setflag) - error->all(FLERR,"Compute bond/local set used with without a variable"); - + error->all(FLERR, "Compute bond/local set used with without a variable"); // set singleflag if need to call bond->single() // set velflag if compute any quantities based on velocities @@ -137,16 +167,20 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : singleflag = 0; velflag = 0; for (int i = 0; i < nvalues; i++) { - if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || - bstyle[i] == FY || bstyle[i] == FZ || bstyle[i] == BN) singleflag = 1; - if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || - bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1; + if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || bstyle[i] == FY || + bstyle[i] == FZ || bstyle[i] == BN) + singleflag = 1; + if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || bstyle[i] == ENGVIB || + bstyle[i] == ENGROT) + velflag = 1; } // initialize output - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -157,13 +191,13 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : ComputeBondLocal::~ComputeBondLocal() { - delete [] bstyle; - delete [] bindex; - for (int i = 0; i < nvar; i++) delete [] vstr[i]; - delete [] vstr; - delete [] vvar; + delete[] bstyle; + delete[] bindex; + for (int i = 0; i < nvar; i++) delete[] vstr[i]; + delete[] vstr; + delete[] vvar; - delete [] dstr; + delete[] dstr; memory->destroy(vlocal); memory->destroy(alocal); @@ -173,8 +207,7 @@ ComputeBondLocal::~ComputeBondLocal() void ComputeBondLocal::init() { - if (force->bond == nullptr) - error->all(FLERR,"No bond style is defined for compute bond/local"); + if (force->bond == nullptr) error->all(FLERR, "No bond style is defined for compute bond/local"); for (int i = 0; i < nvalues; i++) if (bstyle[i] == BN && bindex[i] >= force->bond->single_extra) @@ -183,21 +216,21 @@ void ComputeBondLocal::init() if (nvar) { for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR,"Variable name for compute bond/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute bond/local does not exist"); } if (dstr) { dvar = input->variable->find(dstr); - if (dvar < 0) - error->all(FLERR,"Variable name for compute bond/local does not exist"); + if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist"); } } // set ghostvelflag if need to acquire ghost atom velocities - if (velflag && !comm->ghost_velocity) ghostvelflag = 1; - else ghostvelflag = 0; + if (velflag && !comm->ghost_velocity) + ghostvelflag = 1; + else + ghostvelflag = 0; // do initial memory allocation so that memory_usage() is correct @@ -236,17 +269,17 @@ void ComputeBondLocal::compute_local() int ComputeBondLocal::compute_bonds(int flag) { - int i,m,nb,atom1,atom2,imol,iatom,btype,ivar; + int i, m, nb, atom1, atom2, imol, iatom, btype, ivar; tagint tagprev; - double dx,dy,dz,rsq; - double mass1,mass2,masstotal,invmasstotal; - double xcm[3],vcm[3]; - double delr1[3],delr2[3],delv1[3],delv2[3]; - double r12[3],vpar1,vpar2; - double vvib,vrotsq; - double inertia,omegasq; + double dx, dy, dz, rsq; + double mass1, mass2, masstotal, invmasstotal; + double xcm[3], vcm[3]; + double delr1[3], delr2[3], delv1[3], delv2[3]; + double r12[3], vpar1, vpar2; + double vvib, vrotsq; + double inertia, omegasq; double mvv2e; - double engpot,engtrans,engvib,engrot,fbond; + double engpot, engtrans, engvib, engrot, fbond; double *ptr; double **x = atom->x; @@ -280,7 +313,8 @@ int ComputeBondLocal::compute_bonds(int flag) for (atom1 = 0; atom1 < nlocal; atom1++) { if (!(mask[atom1] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) nb = num_bond[atom1]; + if (molecular == Atom::MOLECULAR) + nb = num_bond[atom1]; else { if (molindex[atom1] < 0) continue; imol = molindex[atom1]; @@ -295,7 +329,7 @@ int ComputeBondLocal::compute_bonds(int flag) } else { tagprev = tag[atom1] - iatom - 1; btype = onemols[imol]->bond_type[iatom][i]; - atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev); + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev); } if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; @@ -310,55 +344,56 @@ int ComputeBondLocal::compute_bonds(int flag) dx = x[atom1][0] - x[atom2][0]; dy = x[atom1][1] - x[atom2][1]; dz = x[atom1][2] - x[atom2][2]; - domain->minimum_image(dx,dy,dz); - rsq = dx*dx + dy*dy + dz*dz; + domain->minimum_image(dx, dy, dz); + rsq = dx * dx + dy * dy + dz * dz; if (btype == 0) { fbond = 0.0; } else { - if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond); - else fbond = engpot = 0.0; + if (singleflag) + engpot = bond->single(btype, rsq, atom1, atom2, fbond); + else + fbond = engpot = 0.0; engvib = engrot = engtrans = omegasq = vvib = 0.0; if (velflag) { if (rmass) { mass1 = rmass[atom1]; mass2 = rmass[atom2]; - } - else { + } else { mass1 = mass[type[atom1]]; mass2 = mass[type[atom2]]; } - masstotal = mass1+mass2; + masstotal = mass1 + mass2; invmasstotal = 1.0 / (masstotal); - xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal; - xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal; - xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal; - vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal; - vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal; - vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal; + xcm[0] = (mass1 * x[atom1][0] + mass2 * x[atom2][0]) * invmasstotal; + xcm[1] = (mass1 * x[atom1][1] + mass2 * x[atom2][1]) * invmasstotal; + xcm[2] = (mass1 * x[atom1][2] + mass2 * x[atom2][2]) * invmasstotal; + vcm[0] = (mass1 * v[atom1][0] + mass2 * v[atom2][0]) * invmasstotal; + vcm[1] = (mass1 * v[atom1][1] + mass2 * v[atom2][1]) * invmasstotal; + vcm[2] = (mass1 * v[atom1][2] + mass2 * v[atom2][2]) * invmasstotal; engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm); // r12 = unit bond vector from atom1 to atom2 - MathExtra::sub3(x[atom2],x[atom1],r12); + MathExtra::sub3(x[atom2], x[atom1], r12); MathExtra::norm3(r12); // delr = vector from COM to each atom // delv = velocity of each atom relative to COM - MathExtra::sub3(x[atom1],xcm,delr1); - MathExtra::sub3(x[atom2],xcm,delr2); - MathExtra::sub3(v[atom1],vcm,delv1); - MathExtra::sub3(v[atom2],vcm,delv2); + MathExtra::sub3(x[atom1], xcm, delr1); + MathExtra::sub3(x[atom2], xcm, delr2); + MathExtra::sub3(v[atom1], vcm, delv1); + MathExtra::sub3(v[atom2], vcm, delv2); // vpar = component of delv parallel to bond vector - vpar1 = MathExtra::dot3(delv1,r12); - vpar2 = MathExtra::dot3(delv2,r12); - engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2); + vpar1 = MathExtra::dot3(delv1, r12); + vpar2 = MathExtra::dot3(delv2, r12); + engvib = 0.5 * (mass1 * vpar1 * vpar1 + mass2 * vpar2 * vpar2); // vvib = relative velocity of 2 atoms along bond direction // vvib < 0 for 2 atoms moving towards each other @@ -369,9 +404,8 @@ int ComputeBondLocal::compute_bonds(int flag) // vrotsq = tangential speed squared of atom1 only // omegasq = omega squared, and is the same for atom1 and atom2 - inertia = mass1*MathExtra::lensq3(delr1) + - mass2*MathExtra::lensq3(delr2); - vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1; + inertia = mass1 * MathExtra::lensq3(delr1) + mass2 * MathExtra::lensq3(delr2); + vrotsq = MathExtra::lensq3(delv1) - vpar1 * vpar1; omegasq = vrotsq / MathExtra::lensq3(delr1); engrot = 0.5 * inertia * omegasq; @@ -384,12 +418,14 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; if (nvar) { ivar = 0; - if (dstr) input->variable->internal_set(dvar,sqrt(rsq)); + if (dstr) input->variable->internal_set(dvar, sqrt(rsq)); } // to make sure dx, dy and dz are always from the lower to the higher id @@ -397,55 +433,55 @@ int ComputeBondLocal::compute_bonds(int flag) for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { - case DIST: - ptr[n] = sqrt(rsq); - break; - case DX: - ptr[n] = dx*directionCorrection; - break; - case DY: - ptr[n] = dy*directionCorrection; - break; - case DZ: - ptr[n] = dz*directionCorrection; - break; - case ENGPOT: - ptr[n] = engpot; - break; - case FORCE: - ptr[n] = sqrt(rsq)*fbond; - break; - case FX: - ptr[n] = dx*fbond; - break; - case FY: - ptr[n] = dy*fbond; - break; - case FZ: - ptr[n] = dz*fbond; - break; - case ENGVIB: - ptr[n] = engvib; - break; - case ENGROT: - ptr[n] = engrot; - break; - case ENGTRANS: - ptr[n] = engtrans; - break; - case OMEGA: - ptr[n] = sqrt(omegasq); - break; - case VELVIB: - ptr[n] = vvib; - break; - case VARIABLE: - ptr[n] = input->variable->compute_equal(vvar[ivar]); - ivar++; - break; - case BN: - ptr[n] = bond->svector[bindex[n]]; - break; + case DIST: + ptr[n] = sqrt(rsq); + break; + case DX: + ptr[n] = dx * directionCorrection; + break; + case DY: + ptr[n] = dy * directionCorrection; + break; + case DZ: + ptr[n] = dz * directionCorrection; + break; + case ENGPOT: + ptr[n] = engpot; + break; + case FORCE: + ptr[n] = sqrt(rsq) * fbond; + break; + case FX: + ptr[n] = dx * fbond; + break; + case FY: + ptr[n] = dy * fbond; + break; + case FZ: + ptr[n] = dz * fbond; + break; + case ENGVIB: + ptr[n] = engvib; + break; + case ENGROT: + ptr[n] = engrot; + break; + case ENGTRANS: + ptr[n] = engtrans; + break; + case OMEGA: + ptr[n] = sqrt(omegasq); + break; + case VELVIB: + ptr[n] = vvib; + break; + case VARIABLE: + ptr[n] = input->variable->compute_equal(vvar[ivar]); + ivar++; + break; + case BN: + ptr[n] = bond->svector[bindex[n]]; + break; } } } @@ -459,10 +495,10 @@ int ComputeBondLocal::compute_bonds(int flag) /* ---------------------------------------------------------------------- */ -int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { - int i,j,m; + int i, j, m; double **v = atom->v; @@ -481,7 +517,7 @@ int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; double **v = atom->v; @@ -504,11 +540,11 @@ void ComputeBondLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"bond/local:vector_local"); + memory->create(vlocal, nmax, "bond/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"bond/local:array_local"); + memory->create(alocal, nmax, nvalues, "bond/local:array_local"); array_local = alocal; } } @@ -519,6 +555,6 @@ void ComputeBondLocal::reallocate(int n) double ComputeBondLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index 894d0e33e4..11e970b8de 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,19 +12,21 @@ ------------------------------------------------------------------------- */ #include "compute_dihedral_local.h" -#include -#include + #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "update.h" #include "domain.h" +#include "error.h" #include "force.h" #include "input.h" -#include "variable.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "molecule.h" +#include "update.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -37,14 +38,13 @@ enum { PHI, VARIABLE }; /* ---------------------------------------------------------------------- */ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr), + vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute dihedral/local command"); if (atom->avec->dihedrals_allow == 0) - error->all(FLERR, - "Compute dihedral/local used when dihedrals are not allowed"); + error->all(FLERR, "Compute dihedral/local used when dihedrals are not allowed"); local_flag = 1; @@ -52,7 +52,7 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = narg - 3; bstyle = new int[nvalues]; - vstr = new char*[nvalues]; + vstr = new char *[nvalues]; vvar = new int[nvalues]; nvalues = 0; @@ -60,13 +60,14 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"phi") == 0) { + if (strcmp(arg[iarg], "phi") == 0) { bstyle[nvalues++] = PHI; - } else if (strncmp(arg[iarg],"v_",2) == 0) { + } else if (strncmp(arg[iarg], "v_", 2) == 0) { bstyle[nvalues++] = VARIABLE; vstr[nvar] = utils::strdup(&arg[iarg][2]); nvar++; - } else break; + } else + break; } // optional args @@ -75,47 +76,45 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : pstr = nullptr; while (iarg < narg) { - if (strcmp(arg[iarg],"set") == 0) { + if (strcmp(arg[iarg], "set") == 0) { setflag = 1; - if (iarg+3 > narg) - error->all(FLERR,"Illegal compute dihedral/local command"); - if (strcmp(arg[iarg+1],"phi") == 0) { - delete [] pstr; - pstr = utils::strdup(arg[iarg+2]); - } else error->all(FLERR,"Illegal compute dihedral/local command"); + if (iarg + 3 > narg) error->all(FLERR, "Illegal compute dihedral/local command"); + if (strcmp(arg[iarg + 1], "phi") == 0) { + delete[] pstr; + pstr = utils::strdup(arg[iarg + 2]); + } else + error->all(FLERR, "Illegal compute dihedral/local command"); iarg += 3; - } else error->all(FLERR,"Illegal compute dihedral/local command"); + } else + error->all(FLERR, "Illegal compute dihedral/local command"); } // error check if (nvar) { - if (!setflag) - error->all(FLERR,"Compute dihedral/local variable requires a set variable"); + if (!setflag) error->all(FLERR, "Compute dihedral/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR, - "Variable name for copute dihedral/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for copute dihedral/local does not exist"); if (!input->variable->equalstyle(vvar[i])) - error->all(FLERR,"Variable for compute dihedral/local is invalid style"); + error->all(FLERR, "Variable for compute dihedral/local is invalid style"); } if (pstr) { pvar = input->variable->find(pstr); - if (pvar < 0) - error->all(FLERR, - "Variable name for compute dihedral/local does not exist"); + if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); if (!input->variable->internalstyle(pvar)) - error->all(FLERR,"Variable for compute dihedral/local is invalid style"); + error->all(FLERR, "Variable for compute dihedral/local is invalid style"); } } else if (setflag) - error->all(FLERR,"Compute dihedral/local set with no variable"); + error->all(FLERR, "Compute dihedral/local set with no variable"); // initialize output - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -126,12 +125,12 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : ComputeDihedralLocal::~ComputeDihedralLocal() { - delete [] bstyle; - for (int i = 0; i < nvar; i++) delete [] vstr[i]; - delete [] vstr; - delete [] vvar; + delete[] bstyle; + for (int i = 0; i < nvar; i++) delete[] vstr[i]; + delete[] vstr; + delete[] vvar; - delete [] pstr; + delete[] pstr; memory->destroy(vlocal); memory->destroy(alocal); @@ -142,21 +141,17 @@ ComputeDihedralLocal::~ComputeDihedralLocal() void ComputeDihedralLocal::init() { if (force->dihedral == nullptr) - error->all(FLERR,"No dihedral style is defined for compute dihedral/local"); + error->all(FLERR, "No dihedral style is defined for compute dihedral/local"); if (nvar) { for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR, - "Variable name for compute dihedral/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); } if (pstr) { pvar = input->variable->find(pstr); - if (pvar < 0) - error->all(FLERR, - "Variable name for compute dihedral/local does not exist"); + if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); } } @@ -191,11 +186,11 @@ void ComputeDihedralLocal::compute_local() int ComputeDihedralLocal::compute_dihedrals(int flag) { - int i,m,nd,atom1,atom2,atom3,atom4,imol,iatom,ivar; + int i, m, nd, atom1, atom2, atom3, atom4, imol, iatom, ivar; tagint tagprev; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv; - double s,c,phi; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, ra2inv, rb2inv, rabinv; + double s, c, phi; double *ptr; double **x = atom->x; @@ -220,7 +215,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) nd = num_dihedral[atom2]; + if (molecular == Atom::MOLECULAR) + nd = num_dihedral[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -237,9 +233,9 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) } 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); + 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; @@ -256,64 +252,66 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) vb1x = x[atom1][0] - x[atom2][0]; vb1y = x[atom1][1] - x[atom2][1]; vb1z = x[atom1][2] - x[atom2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); + domain->minimum_image(vb1x, vb1y, vb1z); vb2x = x[atom3][0] - x[atom2][0]; vb2y = x[atom3][1] - x[atom2][1]; vb2z = x[atom3][2] - x[atom2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); + domain->minimum_image(vb2x, vb2y, vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); + domain->minimum_image(vb2xm, vb2ym, vb2zm); vb3x = x[atom4][0] - x[atom3][0]; vb3y = x[atom4][1] - x[atom3][1]; vb3z = x[atom4][2] - x[atom3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); + domain->minimum_image(vb3x, vb3y, vb3z); - ax = vb1y*vb2zm - vb1z*vb2ym; - ay = vb1z*vb2xm - vb1x*vb2zm; - az = vb1x*vb2ym - vb1y*vb2xm; - bx = vb3y*vb2zm - vb3z*vb2ym; - by = vb3z*vb2xm - vb3x*vb2zm; - bz = vb3x*vb2ym - vb3y*vb2xm; + ax = vb1y * vb2zm - vb1z * vb2ym; + ay = vb1z * vb2xm - vb1x * vb2zm; + az = vb1x * vb2ym - vb1y * vb2xm; + bx = vb3y * vb2zm - vb3z * vb2ym; + by = vb3z * vb2xm - vb3x * vb2zm; + bz = vb3x * vb2ym - vb3y * vb2xm; - rasq = ax*ax + ay*ay + az*az; - rbsq = bx*bx + by*by + bz*bz; - rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm; rg = sqrt(rgsq); ra2inv = rb2inv = 0.0; - if (rasq > 0) ra2inv = 1.0/rasq; - if (rbsq > 0) rb2inv = 1.0/rbsq; - rabinv = sqrt(ra2inv*rb2inv); + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + rabinv = sqrt(ra2inv * rb2inv); - c = (ax*bx + ay*by + az*bz)*rabinv; - s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + c = (ax * bx + ay * by + az * bz) * rabinv; + s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - phi = atan2(s,c); + phi = atan2(s, c); - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; if (nvar) { ivar = 0; - if (pstr) input->variable->internal_set(pvar,phi); + if (pstr) input->variable->internal_set(pvar, phi); } for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { - case PHI: - ptr[n] = 180.0*phi/MY_PI; - break; - case VARIABLE: - ptr[n] = input->variable->compute_equal(vvar[ivar]); - ivar++; - break; + case PHI: + ptr[n] = 180.0 * phi / MY_PI; + break; + case VARIABLE: + ptr[n] = input->variable->compute_equal(vvar[ivar]); + ivar++; + break; } } @@ -334,11 +332,11 @@ void ComputeDihedralLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"dihedral/local:vector_local"); + memory->create(vlocal, nmax, "dihedral/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"dihedral/local:array_local"); + memory->create(alocal, nmax, nvalues, "dihedral/local:array_local"); array_local = alocal; } } @@ -349,6 +347,6 @@ void ComputeDihedralLocal::reallocate(int n) double ComputeDihedralLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_improper_local.cpp b/src/compute_improper_local.cpp index a58f4f4d0d..ff8f069a03 100644 --- a/src/compute_improper_local.cpp +++ b/src/compute_improper_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,36 +12,35 @@ ------------------------------------------------------------------------- */ #include "compute_improper_local.h" -#include -#include + #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "update.h" #include "domain.h" +#include "error.h" #include "force.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "molecule.h" +#include "update.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; static constexpr int DELTA = 10000; - -static constexpr double SMALL = 0.001; +static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute improper/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute improper/local command"); if (atom->avec->impropers_allow == 0) - error->all(FLERR, - "Compute improper/local used when impropers are not allowed"); + error->all(FLERR, "Compute improper/local used when impropers are not allowed"); local_flag = 1; nvalues = narg - 3; @@ -50,12 +48,16 @@ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = 0; for (int iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"chi") == 0) cflag = nvalues++; - else error->all(FLERR,"Invalid keyword in compute improper/local command"); + if (strcmp(arg[iarg], "chi") == 0) + cflag = nvalues++; + else + error->all(FLERR, "Invalid keyword in compute improper/local command"); } - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -75,7 +77,7 @@ ComputeImproperLocal::~ComputeImproperLocal() void ComputeImproperLocal::init() { if (force->improper == nullptr) - error->all(FLERR,"No improper style is defined for compute improper/local"); + error->all(FLERR, "No improper style is defined for compute improper/local"); // do initial memory allocation so that memory_usage() is correct @@ -108,11 +110,11 @@ void ComputeImproperLocal::compute_local() int ComputeImproperLocal::compute_impropers(int flag) { - int i,m,n,ni,atom1,atom2,atom3,atom4,imol,iatom; + int i, m, n, ni, atom1, atom2, atom3, atom4, imol, iatom; tagint tagprev; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z; - double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2; - double s12,c; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z; + double ss1, ss2, ss3, r1, r2, r3, c0, c1, c2, s1, s2; + double s12, c; double *cbuf; double **x = atom->x; @@ -135,8 +137,10 @@ int ComputeImproperLocal::compute_impropers(int flag) if (nvalues == 1) { if (cflag >= 0) cbuf = vlocal; } else { - if (cflag >= 0 && alocal) cbuf = &alocal[0][cflag]; - else cbuf = nullptr; + if (cflag >= 0 && alocal) + cbuf = &alocal[0][cflag]; + else + cbuf = nullptr; } } @@ -144,7 +148,8 @@ int ComputeImproperLocal::compute_impropers(int flag) for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) ni = num_improper[atom2]; + if (molecular == Atom::MOLECULAR) + ni = num_improper[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -161,9 +166,9 @@ int ComputeImproperLocal::compute_impropers(int flag) } else { if (tag[atom2] != onemols[imol]->improper_atom2[atom2][i]) continue; tagprev = tag[atom2] - iatom - 1; - atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i]+tagprev); - atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i]+tagprev); - atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i]+tagprev); + atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i] + tagprev); } if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; @@ -178,21 +183,21 @@ int ComputeImproperLocal::compute_impropers(int flag) vb1x = x[atom1][0] - x[atom2][0]; vb1y = x[atom1][1] - x[atom2][1]; vb1z = x[atom1][2] - x[atom2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); + domain->minimum_image(vb1x, vb1y, vb1z); vb2x = x[atom3][0] - x[atom2][0]; vb2y = x[atom3][1] - x[atom2][1]; vb2z = x[atom3][2] - x[atom2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); + domain->minimum_image(vb2x, vb2y, vb2z); vb3x = x[atom4][0] - x[atom3][0]; vb3y = x[atom4][1] - x[atom3][1]; vb3z = x[atom4][2] - x[atom3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); + domain->minimum_image(vb3x, vb3y, vb3z); - ss1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - ss2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - ss3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + ss1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + ss2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z); + ss3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); r1 = sqrt(ss1); r2 = sqrt(ss2); @@ -202,20 +207,20 @@ int ComputeImproperLocal::compute_impropers(int flag) c1 = (vb1x * vb2x + vb1y * vb2y + vb1z * vb2z) * r1 * r2; c2 = -(vb3x * vb2x + vb3y * vb2y + vb3z * vb2z) * r3 * r2; - s1 = 1.0 - c1*c1; + s1 = 1.0 - c1 * c1; if (s1 < SMALL) s1 = SMALL; s1 = 1.0 / s1; - s2 = 1.0 - c2*c2; + s2 = 1.0 - c2 * c2; if (s2 < SMALL) s2 = SMALL; s2 = 1.0 / s2; - s12 = sqrt(s1*s2); - c = (c1*c2 + c0) * s12; + s12 = sqrt(s1 * s2); + c = (c1 * c2 + c0) * s12; if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - cbuf[n] = 180.0*acos(c)/MY_PI; + cbuf[n] = 180.0 * acos(c) / MY_PI; } n += nvalues; } @@ -237,11 +242,11 @@ void ComputeImproperLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"improper/local:vector_local"); + memory->create(vlocal, nmax, "improper/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"improper/local:array_local"); + memory->create(alocal, nmax, nvalues, "improper/local:array_local"); array_local = alocal; } } @@ -252,6 +257,6 @@ void ComputeImproperLocal::reallocate(int n) double ComputeImproperLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } From ec5dfcd0fb2bfea3ce7c82f4d599d105f1f055ee Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 22 Oct 2024 09:50:30 -0400 Subject: [PATCH 17/24] remove note that was obsoleted 5 years ago --- doc/src/fix_langevin.rst | 6 ------ 1 file changed, 6 deletions(-) diff --git a/doc/src/fix_langevin.rst b/doc/src/fix_langevin.rst index e04805427e..30e4c48270 100644 --- a/doc/src/fix_langevin.rst +++ b/doc/src/fix_langevin.rst @@ -231,12 +231,6 @@ the particles. As described below, this energy can then be printed out or added to the potential energy of the system to monitor energy conservation. -.. note:: - - This accumulated energy does NOT include kinetic energy removed - by the *zero* flag. LAMMPS will print a warning when both options are - active. - The keyword *zero* can be used to eliminate drift due to the thermostat. Because the random forces on different atoms are independent, they do not sum exactly to zero. As a result, this fix From 83db9b8fe62f64267071a3b66006be2d82c3e9a6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 22 Oct 2024 10:25:22 -0400 Subject: [PATCH 18/24] small cleanup and modernization --- src/compute_pair_local.cpp | 4 ++-- src/compute_property_local.cpp | 16 ++++++---------- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 351b499468..57f15264f0 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -38,7 +38,7 @@ enum { TYPE, RADIUS }; ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), pstyle(nullptr), pindex(nullptr), vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR, "Illegal compute pair/local command"); + if (narg < 4) utils::missing_cmd_args(FLERR, "compute pair/local", error); local_flag = 1; nvalues = narg - 3; @@ -97,7 +97,7 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : // error check if (cutstyle == RADIUS && !atom->radius_flag) - error->all(FLERR, "This compute pair/local requires atom attribute radius"); + error->all(FLERR, "Compute pair/local with ID {} requires atom attribute radius", id); // set singleflag if need to call pair->single() diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index 64f3859117..872f4173ec 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -39,7 +39,7 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr), indices(nullptr), pack_choice(nullptr) { - if (narg < 4) error->all(FLERR, "Illegal compute property/local command"); + if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/local", error); local_flag = 1; nvalues = narg - 3; @@ -202,25 +202,23 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg], "cutoff") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Illegal compute property/local command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute property/local cutoff", error); if (strcmp(arg[iarg + 1], "type") == 0) cutstyle = TYPE; else if (strcmp(arg[iarg + 1], "radius") == 0) cutstyle = RADIUS; else - error->all(FLERR, "Illegal compute property/local command"); + error->all(FLERR, "Unknown compute property/local cutoff keyword: {}", arg[iarg + 1]); iarg += 2; } else - error->all(FLERR, "Illegal compute property/local command"); + error->all(FLERR, "Unknown compute property/local keyword: {}", arg[iarg]); } // error check if (atom->molecular == 2 && (kindflag == BOND || kindflag == ANGLE || kindflag == DIHEDRAL || kindflag == IMPROPER)) - error->all(FLERR, - "Compute property/local does not (yet) work " - "with atom_style template"); + error->all(FLERR, "Compute property/local does not (yet) work with atom_style template"); if (kindflag == BOND && atom->avec->bonds_allow == 0) error->all(FLERR, "Compute property/local for property that isn't allocated"); @@ -231,7 +229,7 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : if (kindflag == IMPROPER && atom->avec->impropers_allow == 0) error->all(FLERR, "Compute property/local for property that isn't allocated"); if (cutstyle == RADIUS && !atom->radius_flag) - error->all(FLERR, "Compute property/local requires atom attribute radius"); + error->all(FLERR, "Compute property/local with ID {} requires atom attribute radius", id); nmax = 0; vlocal = nullptr; @@ -255,8 +253,6 @@ void ComputePropertyLocal::init() if (kindflag == NEIGH || kindflag == PAIR) { if (force->pair == nullptr) error->all(FLERR, "No pair style is defined for compute property/local"); - if (force->pair->single_enable == 0) - error->all(FLERR, "Pair style does not support compute property/local"); } // for NEIGH/PAIR need an occasional half neighbor list From 077c77f40267ace14385cba2c425aaba4584f04d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Oct 2024 12:09:35 -0400 Subject: [PATCH 19/24] contributions from bonded interactions is broken when running in parallel see https://matsci.org/t/missing-bond-contributions-from-compute-stress-mop/58455 for details. --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 20 ++++++++++++--- .../compute_stress_mop_profile.cpp | 25 +++++++++++++------ 2 files changed, 34 insertions(+), 11 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index b8d21d2a4f..c631ad8f29 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -242,34 +242,46 @@ void ComputeStressMop::init() // issue an error for unimplemented intramolecular potentials or Kspace. - if (force->bond) bondflag = 1; + if (force->bond) { + bondflag = 1; + if (comm->nprocs > 1) + error->one(FLERR, "compute stress/mop with bonds does not (yet) support MPI parallel runs"); + } + 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 does not account for angle potentials"); + error->one(FLERR, "compute stress/mop does not account for angle potentials"); } else { angleflag = 1; + if (comm->nprocs > 1) + error->one(FLERR, + "compute stress/mop with angles does not (yet) support MPI parallel runs"); } } 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 does not account for dihedral potentials"); + error->one(FLERR, "compute stress/mop does not account for dihedral potentials"); } else { dihedralflag = 1; + if (comm->nprocs > 1) + error->one(FLERR, + "compute stress/mop with dihedrals does not (yet) support MPI parallel runs"); } } if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) - error->all(FLERR, "compute stress/mop does not account for improper potentials"); + error->one(FLERR, "compute stress/mop does not account for improper potentials"); } if (force->kspace) error->warning(FLERR, "compute stress/mop does not account for kspace contributions"); } // need an occasional half neighbor list + neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); } diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index ca2d095fd9..9b2f451d8b 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -229,30 +229,41 @@ void ComputeStressMopProfile::init() if (comm->me == 0) { - // Compute stress/mop/profile only accounts for pair interactions. - // issue an error if any intramolecular potential or Kspace is defined. + // issue an error for unimplemented intramolecular potentials or Kspace. - if (force->bond) bondflag = 1; + if (force->bond) { + bondflag = 1; + if (comm->nprocs > 1) + error->one( + FLERR, + "compute stress/mop/profile with bonds does not (yet) support MPI parallel runs"); + } 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"); + error->one(FLERR, "compute stress/mop/profile does not account for angle potentials"); } else { angleflag = 1; + if (comm->nprocs > 1) + error->one( + FLERR, + "compute stress/mop/profile with angles does not (yet) support MPI parallel runs"); } } - 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"); + error->one(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); } else { dihedralflag = 1; + if (comm->nprocs > 1) + error->one( + FLERR, + "compute stress/mop/profile with dihedrals does not (yet) support MPI parallel runs"); } } - if (force->improper) if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) From 558b1c519744cfb73109f7b9ef0f434e87df9b09 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Oct 2024 12:10:28 -0400 Subject: [PATCH 20/24] apply clang-format, fix other minor formatting issues, use error->one() --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 126 +++++----- .../compute_stress_mop_profile.cpp | 237 +++++++++--------- 2 files changed, 175 insertions(+), 188 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index c631ad8f29..2b46f80a7a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -39,7 +39,7 @@ using namespace LAMMPS_NS; -static constexpr double SMALL = 0.001; +static constexpr double SMALL = 0.001; enum { X, Y, Z }; enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; @@ -63,7 +63,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( } else if (strcmp(arg[3], "z") == 0) { dir = Z; } else - error->all(FLERR, "Illegal compute stress/mop command"); + error->all(FLERR, "Illegal compute stress/mop plane direction {}", arg[3]); // Position of the plane @@ -89,7 +89,7 @@ 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])) { + if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) { pos1 = pos + domain->prd[dir]; } else { pos1 = pos - domain->prd[dir]; @@ -133,27 +133,28 @@ 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++) { + } 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; + error->all(FLERR, "Unknown compute stress/mop keyword {}", arg[iarg]); iarg++; } // Error checks: - // 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"); + // 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; @@ -173,8 +174,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"); + memory->create(dihedral_local, nvalues, "stress/mop:dihedral_local"); + memory->create(dihedral_global, nvalues, "stress/mop:dihedral_global"); size_vector = nvalues; vector_flag = 1; @@ -202,7 +203,6 @@ ComputeStressMop::~ComputeStressMop() void ComputeStressMop::init() { - // Conversion constants nktv2p = force->nktv2p; @@ -224,11 +224,12 @@ void ComputeStressMop::init() dt = update->dt; - // Error check + // Error checks: // Compute stress/mop requires fixed simulation box + if (domain->box_change_size || domain->box_change_shape || domain->deform_flag) - error->all(FLERR, "Compute stress/mop requires a fixed size simulation box"); + error->all(FLERR, "Compute stress/mop requires a fixed simulation box geometry"); // This compute requires a pair style with pair_single method implemented @@ -336,11 +337,11 @@ void ComputeStressMop::compute_vector() //Compute dihedral contribution on separate procs compute_dihedrals(); } else { - for (int i=0; i pos)) || - ((xi[dir] < pos1) && (xj[dir] > pos1))) { + ((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; @@ -881,9 +882,7 @@ void ComputeStressMop::compute_dihedrals() double x_atom_4[3] = {0.0, 0.0, 0.0}; // initialization - for (int i = 0; i < nvalues; i++) { - dihedral_local[i] = 0.0; - } + 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++) { @@ -901,9 +900,9 @@ void ComputeStressMop::compute_dihedrals() 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]); + 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; @@ -955,9 +954,9 @@ void ComputeStressMop::compute_dihedrals() 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 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)); + bool left_cross = ((tau_left >= 0) && (tau_left <= 1)); // no bonds crossing the plane if (!right_cross && !middle_cross && !left_cross) continue; @@ -984,45 +983,45 @@ void ComputeStressMop::compute_dihedrals() 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); + 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; + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; // 1st and 2nd angle - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); + 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); + 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); + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); sc1 = sqrt(sin2); if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; + sc1 = 1.0 / sc1; - sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); sc2 = sqrt(sin2); if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; + sc2 = 1.0 / sc2; s1 = sc1 * sc1; s2 = sc2 * sc2; s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (c0 + c1mag * c2mag) * s12; // error check if (c > 1.0) c = 1.0; @@ -1032,28 +1031,28 @@ void ComputeStressMop::compute_dihedrals() 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); + 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; + 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; + 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; + 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]; @@ -1118,9 +1117,9 @@ void ComputeStressMop::compute_dihedrals() 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; + local_contribution[0] += df[0] / area * nktv2p; + local_contribution[1] += df[1] / area * nktv2p; + local_contribution[2] += df[2] / area * nktv2p; } } @@ -1128,11 +1127,10 @@ void ComputeStressMop::compute_dihedrals() 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]; + dihedral_local[m] = local_contribution[0]; + dihedral_local[m + 1] = local_contribution[1]; + dihedral_local[m + 2] = local_contribution[2]; } m += 3; } - } diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 9b2f451d8b..154cbbcfde 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -39,7 +39,7 @@ using namespace LAMMPS_NS; -static constexpr double SMALL = 0.001; +static constexpr double SMALL = 0.001; enum { X, Y, Z }; enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; @@ -64,7 +64,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a } else if (strcmp(arg[3], "z") == 0) { dir = Z; } else - error->all(FLERR, "Illegal compute stress/mop/profile command"); + error->all(FLERR, "Illegal compute stress/mop/profile plane direction {}", arg[3]); // bin parameters @@ -118,13 +118,13 @@ 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++) { + } 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; + error->all(FLERR, "Unknown compute stress/mop/profile keyword {}", arg[iarg]); iarg++; } @@ -134,14 +134,14 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a // 3D only if (domain->dimension == 2 && dir == Z) - error->all(FLERR, "Compute stress/mop/profile incompatible with Z in 2d system"); + error->all(FLERR, "Compute stress/mop/profile is incompatible with Z in 2d system"); // orthogonal simulation box if (domain->triclinic != 0) - error->all(FLERR, "Compute stress/mop/profile incompatible with triclinic simulation box"); + error->all(FLERR, "Compute stress/mop/profile is incompatible with triclinic simulation box"); - // initialize some variables + // Initialize some variables nbins = 0; coord = coordp = nullptr; @@ -192,12 +192,13 @@ ComputeStressMopProfile::~ComputeStressMopProfile() void ComputeStressMopProfile::init() { - // conversion constants + // Conversion constants nktv2p = force->nktv2p; ftm2v = force->ftm2v; - // plane area + // Plane area + if (domain->dimension == 3) { area = 1; int i; @@ -208,7 +209,7 @@ void ComputeStressMopProfile::init() area = (dir == X) ? domain->prd[1] : domain->prd[0]; } - // timestep Value + // Timestep Value dt = update->dt; @@ -217,13 +218,14 @@ void ComputeStressMopProfile::init() // Compute stress/mop/profile requires fixed simulation box if (domain->box_change_size || domain->box_change_shape || domain->deform_flag) - error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box"); + error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box geometry"); // This compute requires a pair style with pair_single method implemented if (!force->pair) error->all(FLERR, "No pair style is defined for compute stress/mop/profile"); if (force->pair->single_enable == 0) - error->all(FLERR, "Pair style does not support compute stress/mop/profile"); + error->all(FLERR, "Pair style {} does not support compute stress/mop/profile", + force->pair_style); // Errors @@ -264,10 +266,11 @@ void ComputeStressMopProfile::init() "compute stress/mop/profile with dihedrals does not (yet) support MPI parallel runs"); } } - if (force->improper) + if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) - error->all(FLERR, "compute stress/mop/profile does not account for improper potentials"); + error->one(FLERR, "compute stress/mop/profile does not account for improper potentials"); + } if (force->kspace) error->warning(FLERR, "compute stress/mop/profile does not account for kspace contributions"); } @@ -319,28 +322,26 @@ void ComputeStressMopProfile::compute_array() compute_angles(); } else { for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - angle_local[m][i] = 0.0; - } + 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); + 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 i = 0; i < nvalues; i++) { dihedral_local[m][i] = 0.0; } } } // sum dihedral contribution over all procs - MPI_Allreduce(&dihedral_local[0][0],&dihedral_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + 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]; @@ -348,7 +349,8 @@ 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] + dihedral_global[ibin][m]; + array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] + + dihedral_global[ibin][m]; m++; } } @@ -737,15 +739,12 @@ void ComputeStressMopProfile::compute_angles() // initialization for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - angle_local[m][i] = 0.0; - } + 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; @@ -776,7 +775,7 @@ void ComputeStressMopProfile::compute_angles() if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; if (atype <= 0) continue; - for (int ibin = 0; ibin=0) && (tau_right <= 1)); - bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + 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); + 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; 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) - { + 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; + 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) - { + 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; + 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) - { + 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; + 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; + 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; + 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; } } } @@ -874,8 +873,8 @@ void ComputeStressMopProfile::compute_angles() 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]; + angle_local[ibin][m + 1] = local_contribution[ibin][1]; + angle_local[ibin][m + 2] = local_contribution[ibin][2]; } } m += 3; @@ -928,9 +927,7 @@ void ComputeStressMopProfile::compute_dihedrals() // initialization for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - dihedral_local[m][i] = 0.0; - } + 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; @@ -951,9 +948,9 @@ void ComputeStressMopProfile::compute_dihedrals() 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]); + 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; @@ -966,7 +963,7 @@ void ComputeStressMopProfile::compute_dihedrals() if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; - for (int ibin = 0; ibin=0) && (tau_right <= 1)); + 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)); + bool left_cross = ((tau_left >= 0) && (tau_left <= 1)); // no bonds crossing the plane if (!right_cross && !middle_cross && !left_cross) continue; @@ -1037,45 +1034,45 @@ void ComputeStressMopProfile::compute_dihedrals() 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); + 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; + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; // 1st and 2nd angle - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); + 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); + 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); + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); sc1 = sqrt(sin2); if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; + sc1 = 1.0 / sc1; - sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); sc2 = sqrt(sin2); if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; + sc2 = 1.0 / sc2; s1 = sc1 * sc1; s2 = sc2 * sc2; s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (c0 + c1mag * c2mag) * s12; // error check if (c > 1.0) c = 1.0; @@ -1085,36 +1082,35 @@ void ComputeStressMopProfile::compute_dihedrals() 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); + 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; + 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; + 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; + 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) - { + 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]; @@ -1122,8 +1118,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only middle bond crossing the plane - if (!right_cross && middle_cross && !left_cross) - { + 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]); @@ -1131,8 +1126,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only left bond crossing the plane - if (!right_cross && !middle_cross && left_cross) - { + 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]; @@ -1140,8 +1134,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only right & middle bonds crossing the plane - if (right_cross && middle_cross && !left_cross) - { + 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]; @@ -1149,8 +1142,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only right & left bonds crossing the plane - if (right_cross && !middle_cross && left_cross) - { + 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]); @@ -1158,8 +1150,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only middle & left bonds crossing the plane - if (!right_cross && middle_cross && left_cross) - { + 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]; @@ -1167,17 +1158,16 @@ void ComputeStressMopProfile::compute_dihedrals() } // all three bonds crossing the plane - if (right_cross && middle_cross && left_cross) - { + 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; + local_contribution[ibin][0] += df[0] / area * nktv2p; + local_contribution[ibin][1] += df[1] / area * nktv2p; + local_contribution[ibin][2] += df[2] / area * nktv2p; } } } @@ -1188,13 +1178,12 @@ void ComputeStressMopProfile::compute_dihedrals() 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]; + dihedral_local[ibin][m + 1] = local_contribution[ibin][1]; + dihedral_local[ibin][m + 2] = local_contribution[ibin][2]; } } m += 3; } - } /* ---------------------------------------------------------------------- @@ -1214,11 +1203,11 @@ void ComputeStressMopProfile::setup_bins() if ((origin > domain->boxhi[dir]) || (origin < domain->boxlo[dir])) 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((origin - boxlo[dir]) * invdelta); + lo = origin - n * delta; - n = static_cast ((boxhi[dir] - origin) * invdelta); - hi = origin + n*delta; + n = static_cast((boxhi[dir] - origin) * invdelta); + hi = origin + n * delta; offset = lo; nbins = static_cast((hi - lo) * invdelta + 1.5); @@ -1232,8 +1221,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(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 b810053d0211b0099fb4353f695f95e39487251e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Oct 2024 12:35:50 -0400 Subject: [PATCH 21/24] document restrictions for bonded interactions --- doc/src/compute_stress_mop.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index b4779b8bf3..31c81b5df7 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -129,6 +129,9 @@ package ` doc page on for more info. The method is implemented for orthogonal simulation boxes whose size does not change in time, and axis-aligned planes. +Contributions from bonds, angles, and dihedrals are not compatible +with MPI parallel runs. + 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 From b682f52d111b4f5fbcf504620302e0d4b9d10ddd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Oct 2024 21:07:39 -0400 Subject: [PATCH 22/24] fix typos --- doc/src/Howto_bioFF.rst | 8 ++++---- src/compute_angle_local.cpp | 2 +- src/compute_bond_local.cpp | 2 +- src/compute_dihedral_local.cpp | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index d35913eb0d..92dd45b9b6 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -248,11 +248,11 @@ exist for united atom representations (OPLS-UA) and for all-atom representations (OPLS-AA). This force field is based on atom types mapped to specific functional -groups in organic and biological molecules. Each atom includes an +groups in organic and biological molecules. Each atom includes a static, partial atomic charge reflecting the oxidation state of the -element derived from its bonded neighbors :ref:`(Jorgensen) ` -and computed based on increments determined by the atom type of the -atoms bond to it. +element derived from its bonded neighbors :ref:`(Jorgensen) +` and computed based on increments determined by the +atom type of the atoms bond to it. The interaction styles listed below compute force field formulas that are fully or in part consistent with the OPLS style force fields. See diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index 569c30772e..b4a9334416 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -101,7 +101,7 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : if (!setflag) error->all(FLERR, "Compute angle/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) error->all(FLERR, "Variable name for copute angle/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); if (!input->variable->equalstyle(vvar[i])) error->all(FLERR, "Variable for compute angle/local is invalid style"); } diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index bd59d76b9c..e9632d254f 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -147,7 +147,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); if (vvar[i] < 0) - error->all(FLERR, "Variable name {} for copute bond/local does not exist", vstr[i]); + error->all(FLERR, "Variable name {} for compute bond/local does not exist", vstr[i]); if (!input->variable->equalstyle(vvar[i])) error->all(FLERR, "Variable {} for compute bond/local is invalid style", vstr[i]); } diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index 11e970b8de..bd1e389e85 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -95,7 +95,7 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : if (!setflag) error->all(FLERR, "Compute dihedral/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) error->all(FLERR, "Variable name for copute dihedral/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); if (!input->variable->equalstyle(vvar[i])) error->all(FLERR, "Variable for compute dihedral/local is invalid style"); } From fdd2fc4f5db86afe98c0abb75a707af2d4ad70c8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Oct 2024 21:08:07 -0400 Subject: [PATCH 23/24] move error check to Compute::init() so they cannot be bypassed by input commands --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 22 ++++++++--------- .../compute_stress_mop_profile.cpp | 24 +++++++++---------- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 2b46f80a7a..7c66c1fe50 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -144,17 +144,6 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( iarg++; } - // Error checks: - - // 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"); - - // 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; @@ -203,6 +192,17 @@ ComputeStressMop::~ComputeStressMop() void ComputeStressMop::init() { + // Error checks: + + // 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"); + + // orthogonal simulation box + + if (domain->triclinic != 0) + error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box"); + // Conversion constants nktv2p = force->nktv2p; diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 154cbbcfde..c8087b60a9 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -129,18 +129,6 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a iarg++; } - // check domain related errors - - // 3D only - - if (domain->dimension == 2 && dir == Z) - error->all(FLERR, "Compute stress/mop/profile is incompatible with Z in 2d system"); - - // orthogonal simulation box - - if (domain->triclinic != 0) - error->all(FLERR, "Compute stress/mop/profile is incompatible with triclinic simulation box"); - // Initialize some variables nbins = 0; @@ -192,6 +180,18 @@ ComputeStressMopProfile::~ComputeStressMopProfile() void ComputeStressMopProfile::init() { + // check domain related errors + + // 3D only + + if (domain->dimension == 2 && dir == Z) + error->all(FLERR, "Compute stress/mop/profile is incompatible with Z in 2d system"); + + // check for orthogonal simulation box + + if (domain->triclinic != 0) + error->all(FLERR, "Compute stress/mop/profile is incompatible with triclinic simulation box"); + // Conversion constants nktv2p = force->nktv2p; From 73c049459b4dc6dd6f46d92f5d60693cf5983c2c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Oct 2024 21:10:34 -0400 Subject: [PATCH 24/24] spelling --- doc/src/fix_atom_swap.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_atom_swap.rst b/doc/src/fix_atom_swap.rst index 1a833f69b1..fd4a2f7245 100644 --- a/doc/src/fix_atom_swap.rst +++ b/doc/src/fix_atom_swap.rst @@ -181,8 +181,8 @@ built with that package. See the :doc:`Build package ` doc page for more info. When this fix is used with a :doc:`hybrid pair style ` -system, only swaps between atom types of the same substyle (or -combination of substyles) are permitted. +system, only swaps between atom types of the same sub-style (or +combination of sub-styles) are permitted. This fix cannot be used with systems that do not have per-type masses (e.g. atom style sphere) since the implemented algorithm pre-computes