From e3dd2908d9ac1c1e944b8c648a00ce6acb0c9a17 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 6 Jan 2022 19:42:45 -0500 Subject: [PATCH 01/75] Step version strings for the next patch release --- doc/lammps.1 | 2 +- src/version.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/lammps.1 b/doc/lammps.1 index a374e62604..58086b1fae 100644 --- a/doc/lammps.1 +++ b/doc/lammps.1 @@ -1,4 +1,4 @@ -.TH LAMMPS "1" "14 December 2021" "2021-12-14" +.TH LAMMPS "1" "7 January 2022" "2022-1-7" .SH NAME .B LAMMPS \- Molecular Dynamics Simulator. diff --git a/src/version.h b/src/version.h index 6bb32291be..c33ca6ee15 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "14 Dec 2021" +#define LAMMPS_VERSION "7 Jan 2022" From 240db21054ff8e5445d709caa825d76ee6e3ef57 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 11 Jan 2022 11:46:08 -0500 Subject: [PATCH 02/75] silence possible warnings about missing files on "make clean-all" --- src/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile b/src/Makefile index 2d651c4986..d23058447a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -439,7 +439,7 @@ clean: clean-all: rm -rf Obj_* - rm style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h + rm -f style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h clean-%: @if [ $@ = "clean-serial" ]; \ From 698256f4fea42bcbbaa350f6e64163a32d3a1adf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 12 Jan 2022 08:17:54 -0500 Subject: [PATCH 03/75] update fedora singularity image to Fedora 35 --- tools/singularity/{fedora34_mingw.def => fedora35_mingw.def} | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) rename tools/singularity/{fedora34_mingw.def => fedora35_mingw.def} (99%) diff --git a/tools/singularity/fedora34_mingw.def b/tools/singularity/fedora35_mingw.def similarity index 99% rename from tools/singularity/fedora34_mingw.def rename to tools/singularity/fedora35_mingw.def index 8a1f108aec..c2108a378c 100644 --- a/tools/singularity/fedora34_mingw.def +++ b/tools/singularity/fedora35_mingw.def @@ -1,5 +1,5 @@ BootStrap: docker -From: fedora:34 +From: fedora:35 %post dnf -y update @@ -94,7 +94,7 @@ EOF CUSTOM_PROMPT_ENV=/.singularity.d/env/99-zz_custom_prompt.sh cat >$CUSTOM_PROMPT_ENV < Date: Fri, 14 Jan 2022 13:02:15 -0500 Subject: [PATCH 04/75] check return values for errors --- src/MANYBODY/pair_eam_cd.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp index 46d439b879..dc1ec360c6 100644 --- a/src/MANYBODY/pair_eam_cd.cpp +++ b/src/MANYBODY/pair_eam_cd.cpp @@ -499,10 +499,13 @@ void PairEAMCD::read_h_coeff(char *filename) // Seek to end of file, read last part into a buffer and // then skip over lines in buffer until reaching the end. - platform::fseek(fptr, platform::END_OF_FILE); - platform::fseek(fptr, platform::ftell(fptr) - MAXLINE); + if ( (platform::fseek(fptr, platform::END_OF_FILE) < 0) + || (platform::fseek(fptr, platform::ftell(fptr) - MAXLINE) < 0)) + error->one(FLERR,"Failure to seek to end-of-file for reading h(x) coeffs: {}", + utils::getsyserror()); + char *buf = new char[MAXLINE+1]; - fread(buf, 1, MAXLINE, fptr); + utils::sfread(FLERR, buf, 1, MAXLINE, fptr, filename, error); buf[MAXLINE] = '\0'; // must 0-terminate buffer for string processing Tokenizer lines(buf, "\n"); delete[] buf; From 2213eb8d3f090e88d06a977423e8d240e7201d35 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jan 2022 14:45:12 -0500 Subject: [PATCH 05/75] use enum with symbolic constants instead of numbers --- src/output.cpp | 49 +++++++++++++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/src/output.cpp b/src/output.cpp index a8fda4173d..4156c6fc5e 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -41,6 +41,8 @@ using namespace LAMMPS_NS; #define DELTA 1 #define EPSDT 1.0e-6 +enum {SETUP, WRITE, RESET_DT}; + /* ---------------------------------------------------------------------- initialize all output ------------------------------------------------------------------------- */ @@ -232,7 +234,7 @@ void Output::setup(int memflag) // only do this if dump written or dump has not been written yet if (writeflag || last_dump[idump] < 0) - calculate_next_dump(0,idump,ntimestep); + calculate_next_dump(SETUP,idump,ntimestep); // if dump not written now, use addstep_compute_all() // since don't know what computes the dump will invoke @@ -361,7 +363,7 @@ void Output::setup(int memflag) dump[idump]->write(); last_dump[idump] = ntimestep; - calculate_next_dump(1,idump,ntimestep); + calculate_next_dump(WRITE,idump,ntimestep); if (mode_dump[idump] == 0 && (dump[idump]->clearstep || var_dump[idump])) @@ -468,9 +470,9 @@ void Output::setup(int memflag) for timestep mode, set next_dump for simulation time mode, set next_time_dump and next_dump which flag depends on caller - 0 = from setup() at start of run - 1 = from write() during run each time a dump file is written - 2 = from reset_dt() called from fix dt/reset when it changes timestep size + SETUP = from setup() at start of run + WRITE = from write() during run each time a dump file is written + RESET_DT = from reset_dt() called from fix dt/reset when it changes timestep size ------------------------------------------------------------------------- */ void Output::calculate_next_dump(int which, int idump, bigint ntimestep) @@ -482,13 +484,13 @@ void Output::setup(int memflag) if (every_dump[idump]) { - // which = 0: nextdump = next multiple of every_dump - // which = 1: increment nextdump by every_dump + // which = SETUP: nextdump = next multiple of every_dump + // which = WRITE: increment nextdump by every_dump - if (which == 0) + if (which == SETUP) next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; - else if (which == 1) + else if (which == WRITE) next_dump[idump] += every_dump[idump]; } else { @@ -510,17 +512,28 @@ void Output::setup(int memflag) if (every_time_dump[idump] > 0.0) { - // which = 0: nexttime = next multiple of every_time_dump - // which = 1: increment nexttime by every_time_dump - // which = 2: no change to previous nexttime (only timestep has changed) + // which = SETUP: nexttime = next multiple of every_time_dump + // which = WRITE: increment nexttime by every_time_dump + // which = RESET_DT: no change to previous nexttime (only timestep has changed) - if (which == 0) + switch (which) { + case SETUP: nexttime = static_cast (tcurrent/every_time_dump[idump]) * every_time_dump[idump] + every_time_dump[idump]; - else if (which == 1) + break; + + case WRITE: nexttime = next_time_dump[idump] + every_time_dump[idump]; - else if (which == 2) + break; + + case RESET_DT: nexttime = next_time_dump[idump]; + break; + + default: + nexttime = 0; + error->all(FLERR,"Unexpected argument to calculate_next_dump"); + } nextdump = ntimestep + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / @@ -541,10 +554,10 @@ void Output::setup(int memflag) } else { - // do not re-evaulate variable for which = 2, leave nexttime as-is + // do not re-evaulate variable for which = RESET_DT, leave nexttime as-is // unless next_time_dump < 0.0, which means variable never yet evaluated - if (which < 2 || next_time_dump[idump] < 0.0) { + if (which < RESET_DT || next_time_dump[idump] < 0.0) { nexttime = input->variable->compute_equal(ivar_dump[idump]); } else nexttime = next_time_dump[idump]; @@ -704,7 +717,7 @@ void Output::reset_dt() // since timestep change affects next step if (next_dump[idump] != ntimestep) - calculate_next_dump(2,idump,update->ntimestep); + calculate_next_dump(RESET_DT,idump,update->ntimestep); if (dump[idump]->clearstep || var_dump[idump]) next_time_dump_any = MIN(next_time_dump_any,next_dump[idump]); From dcb1ddb28272e87d76058c5bc2b7c2d9d90ca48c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jan 2022 14:55:13 -0500 Subject: [PATCH 06/75] remove redundant code --- src/fix_move.cpp | 63 ++++++++++++++++++------------------------------ 1 file changed, 24 insertions(+), 39 deletions(-) diff --git a/src/fix_move.cpp b/src/fix_move.cpp index 72fd9b75d2..f7bc4d3640 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -522,7 +522,7 @@ void FixMove::initial_integrate(int /*vflag*/) } } - // for wiggle: X = X0 + A sin(w*dt) + // for wiggle: X = X0 + A sin(w*dt) } else if (mstyle == WIGGLE) { double arg = omega_rotate * delta; @@ -578,19 +578,19 @@ void FixMove::initial_integrate(int /*vflag*/) } } - // for rotate by right-hand rule around omega: - // P = point = vector = point of rotation - // R = vector = axis of rotation - // w = omega of rotation (from period) - // X0 = xoriginal = initial coord of atom - // R0 = runit = unit vector for R - // D = X0 - P = vector from P to X0 - // C = (D dot R0) R0 = projection of atom coord onto R line - // A = D - C = vector from R line to X0 - // B = R0 cross A = vector perp to A in plane of rotation - // A,B define plane of circular rotation around R line - // X = P + C + A cos(w*dt) + B sin(w*dt) - // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) + // for rotate by right-hand rule around omega: + // P = point = vector = point of rotation + // R = vector = axis of rotation + // w = omega of rotation (from period) + // X0 = xoriginal = initial coord of atom + // R0 = runit = unit vector for R + // D = X0 - P = vector from P to X0 + // C = (D dot R0) R0 = projection of atom coord onto R line + // A = D - C = vector from R line to X0 + // B = R0 cross A = vector perp to A in plane of rotation + // A,B define plane of circular rotation around R line + // X = P + C + A cos(w*dt) + B sin(w*dt) + // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) } else if (mstyle == ROTATE) { double arg = omega_rotate * delta; @@ -707,10 +707,10 @@ void FixMove::initial_integrate(int /*vflag*/) } } - // for variable: compute x,v from variables - // NOTE: also allow for changes to extra attributes? - // omega, angmom, theta, quat - // only necessary if prescribed motion involves rotation + // for variable: compute x,v from variables + // NOTE: also allow for changes to extra attributes? + // omega, angmom, theta, quat + // only necessary if prescribed motion involves rotation } else if (mstyle == VARIABLE) { @@ -778,21 +778,16 @@ void FixMove::initial_integrate(int /*vflag*/) } else if (vxvarstr) { if (vxvarstyle == EQUAL) v[i][0] = vx; else v[i][0] = velocity[i][0]; - if (rmass) { - x[i][0] += dtv * v[i][0]; - } else { - x[i][0] += dtv * v[i][0]; - } + x[i][0] += dtv * v[i][0]; } else { if (rmass) { dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; - x[i][0] += dtv * v[i][0]; } else { dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; - x[i][0] += dtv * v[i][0]; } + x[i][0] += dtv * v[i][0]; } if (yvarstr && vyvarstr) { @@ -806,21 +801,16 @@ void FixMove::initial_integrate(int /*vflag*/) } else if (vyvarstr) { if (vyvarstyle == EQUAL) v[i][1] = vy; else v[i][1] = velocity[i][1]; - if (rmass) { - x[i][1] += dtv * v[i][1]; - } else { - x[i][1] += dtv * v[i][1]; - } + x[i][1] += dtv * v[i][1]; } else { if (rmass) { dtfm = dtf / rmass[i]; v[i][1] += dtfm * f[i][1]; - x[i][1] += dtv * v[i][1]; } else { dtfm = dtf / mass[type[i]]; v[i][1] += dtfm * f[i][1]; - x[i][1] += dtv * v[i][1]; } + x[i][1] += dtv * v[i][1]; } if (zvarstr && vzvarstr) { @@ -834,21 +824,16 @@ void FixMove::initial_integrate(int /*vflag*/) } else if (vzvarstr) { if (vzvarstyle == EQUAL) v[i][2] = vz; else v[i][2] = velocity[i][2]; - if (rmass) { - x[i][2] += dtv * v[i][2]; - } else { - x[i][2] += dtv * v[i][2]; - } + x[i][2] += dtv * v[i][2]; } else { if (rmass) { dtfm = dtf / rmass[i]; v[i][2] += dtfm * f[i][2]; - x[i][2] += dtv * v[i][2]; } else { dtfm = dtf / mass[type[i]]; v[i][2] += dtfm * f[i][2]; - x[i][2] += dtv * v[i][2]; } + x[i][2] += dtv * v[i][2]; } domain->remap_near(x[i],xold); From ed702b930928584017f70a5d83a4b32304e8e1d5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jan 2022 15:58:57 -0500 Subject: [PATCH 07/75] don't allow exceptions to "escape" a destructor --- src/REAXFF/fix_reaxff_species.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 4f44cc7c64..bed7fabb20 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -246,8 +246,10 @@ FixReaxFFSpecies::~FixReaxFFSpecies() if (posflag && multipos_opened) fclose(pos); } - modify->delete_compute(fmt::format("SPECATOM_{}",id)); - modify->delete_fix(fmt::format("SPECBOND_{}",id)); + try { + modify->delete_compute(fmt::format("SPECATOM_{}",id)); + modify->delete_fix(fmt::format("SPECBOND_{}",id)); + } catch (std::exception &) {} } /* ---------------------------------------------------------------------- */ From 11cc8a6a5928548a186f7b5536d00ff186f5ded0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jan 2022 15:59:05 -0500 Subject: [PATCH 08/75] whitespace --- src/output.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/output.cpp b/src/output.cpp index 4156c6fc5e..6c14900ade 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -475,7 +475,7 @@ void Output::setup(int memflag) RESET_DT = from reset_dt() called from fix dt/reset when it changes timestep size ------------------------------------------------------------------------- */ - void Output::calculate_next_dump(int which, int idump, bigint ntimestep) +void Output::calculate_next_dump(int which, int idump, bigint ntimestep) { // dump mode is by timestep // just set next_dump From e271a54802014cefd7e8d71de3fa9e125540e84a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jan 2022 20:25:00 -0500 Subject: [PATCH 09/75] re-enable using rerun after the changes from PR #3052 --- src/output.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/output.cpp b/src/output.cpp index a8fda4173d..b2b8011842 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -619,9 +619,8 @@ int Output::check_time_dumps(bigint ntimestep) { next_dump_any = MAXBIGINT; for (int idump = 0; idump < ndump; idump++) - if (last_dump[idump] >= 0) - error->all(FLERR, - "Cannot reset timestep with active dump - must undump first"); + if ((last_dump[idump] >= 0) && !update->whichflag) + error->all(FLERR, "Cannot reset timestep with active dump - must undump first"); if (restart_flag_single) { if (restart_every_single) { From 1e7969c7b9bf3977b08483fb56006c7b607708d3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 15 Jan 2022 19:34:47 -0500 Subject: [PATCH 10/75] add new pair styles harmonic/cut and harmonic/cut/omp --- doc/src/Commands_pair.rst | 1 + doc/src/pair_harmonic_cut.rst | 90 +++++ doc/src/pair_style.rst | 1 + src/.gitignore | 2 + src/EXTRA-PAIR/pair_harmonic_cut.cpp | 314 ++++++++++++++++++ src/EXTRA-PAIR/pair_harmonic_cut.h | 72 ++++ src/OPENMP/pair_harmonic_cut_omp.cpp | 152 +++++++++ src/OPENMP/pair_harmonic_cut_omp.h | 48 +++ .../tests/mol-pair-harmonic_cut.yaml | 94 ++++++ 9 files changed, 774 insertions(+) create mode 100644 doc/src/pair_harmonic_cut.rst create mode 100644 src/EXTRA-PAIR/pair_harmonic_cut.cpp create mode 100644 src/EXTRA-PAIR/pair_harmonic_cut.h create mode 100644 src/OPENMP/pair_harmonic_cut_omp.cpp create mode 100644 src/OPENMP/pair_harmonic_cut_omp.h create mode 100644 unittest/force-styles/tests/mol-pair-harmonic_cut.yaml diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 9ac4fc851c..a226d73b36 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -119,6 +119,7 @@ OPT. * :doc:`granular ` * :doc:`gw ` * :doc:`gw/zbl ` + * :doc:`harmonic/cut (o) ` * :doc:`hbond/dreiding/lj (o) ` * :doc:`hbond/dreiding/morse (o) ` * :doc:`hdnnp ` diff --git a/doc/src/pair_harmonic_cut.rst b/doc/src/pair_harmonic_cut.rst new file mode 100644 index 0000000000..877d2d97ff --- /dev/null +++ b/doc/src/pair_harmonic_cut.rst @@ -0,0 +1,90 @@ +.. index:: pair_style harmonic/cut +.. index:: pair_style harmonic/cut/omp + +pair_style harmonic/cut command +=============================== + +Accelerator Variants: *harmonic/cut/omp* + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style style + +* style = *harmonic/cut* + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style harmonic/cut 2.5 + pair_coeff * * 0.2 2.0 + pair_coeff 1 1 0.5 2.5 + +Description +""""""""""" + +Style *harmonic/cut* computes pairwise repulsive-only harmonic interactions with the formula + +.. math:: + + E = k (r_c - r)^2 \qquad r < r_c + +:math:`r_c` is the cutoff. + +The following coefficients must be defined for each pair of atoms +types via the :doc:`pair_coeff ` command as in the examples +above, or in the data file or restart files read by the +:doc:`read_data ` or :doc:`read_restart ` +commands: + +* :math:`k` (energy units) +* :math:`r_c` (distance units) + +---------- + +.. include:: accel_styles.rst + +---------- + +Mixing, shift, table, tail correction, restart, rRESPA info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +For atom type pairs I,J and I != J, the :math:`k` and :math:`r_c` +coefficients can be mixed. The default mix value is *geometric*. +See the "pair_modify" command for details. + +Since the potential is zero at and beyond the cutoff parameter by +construction, there is no need to support support the :doc:`pair_modify +` shift or tail options for the energy and pressure of the +pair interaction. + +These pair styles write their information to :doc:`binary restart files `, +so pair_style and pair_coeff commands do not need to be specified in an input script +that reads a restart file. + +These pair styles can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. They do not support the +*inner*, *middle*, *outer* keywords. + +---------- + +Restrictions +"""""""""""" + +The *harmonic/cut* pair style is only enabled if LAMMPS was +built with the EXTRA-PAIR package. +See the :doc:`Build package ` page for more info. + +Related commands +"""""""""""""""" + +:doc:`pair_coeff ` + +Default +""""""" + +none diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 4bb3c90a8d..f657e29aa3 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -183,6 +183,7 @@ accelerated styles exist. * :doc:`gran/hooke/history ` - granular potential without history effects * :doc:`gw ` - Gao-Weber potential * :doc:`gw/zbl ` - Gao-Weber potential with a repulsive ZBL core +* :doc:`harmonic/cut ` - repulsive-only harmonic potential * :doc:`hbond/dreiding/lj ` - DREIDING hydrogen bonding LJ potential * :doc:`hbond/dreiding/morse ` - DREIDING hydrogen bonding Morse potential * :doc:`hdnnp ` - High-dimensional neural network potential diff --git a/src/.gitignore b/src/.gitignore index 695c9a19af..a9eef4b371 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1489,6 +1489,8 @@ /pair_dpd_fdt.h /pair_dpd_fdt_energy.cpp /pair_dpd_fdt_energy.h +/pair_harmonic_cut.cpp +/pair_harmonic_cut.h /pair_lennard_mdf.cpp /pair_lennard_mdf.h /pair_lj_cut_coul_long_cs.cpp diff --git a/src/EXTRA-PAIR/pair_harmonic_cut.cpp b/src/EXTRA-PAIR/pair_harmonic_cut.cpp new file mode 100644 index 0000000000..bd0e48efdc --- /dev/null +++ b/src/EXTRA-PAIR/pair_harmonic_cut.cpp @@ -0,0 +1,314 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_harmonic_cut.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairHarmonicCut::PairHarmonicCut(LAMMPS *lmp) : Pair(lmp) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairHarmonicCut::~PairHarmonicCut() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(cut); + memory->destroy(cutsq); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairHarmonicCut::compute(int eflag, int vflag) +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, fxtmp, fytmp, fztmp; + double delx, dely, delz, rsq, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; + + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + fxtmp = fytmp = fztmp = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + const double r = sqrt(rsq); + const double delta = cut[itype][jtype] - r; + const double philj = factor_lj * delta * k[itype][jtype]; + const double fpair = delta * philj / r; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + } + + if (evflag) ev_tally(i, j, nlocal, newton_pair, philj, 0.0, fpair, delx, dely, delz); + } + } + f[i][0] += fxtmp; + f[i][1] += fytmp; + f[i][2] += fztmp; + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairHarmonicCut::allocate() +{ + allocated = 1; + int n = atom->ntypes + 1; + + memory->create(setflag, n, n, "pair:setflag"); + for (int i = 1; i < n; i++) + for (int j = i; j < n; j++) setflag[i][j] = 0; + + memory->create(k, n, n, "pair:k"); + memory->create(cut, n, n, "pair:cut"); + memory->create(cutsq, n, n, "pair:cutsq"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairHarmonicCut::settings(int narg, char **arg) +{ + if (narg > 0) error->all(FLERR, "Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHarmonicCut::coeff(int narg, char **arg) +{ + if (narg != 4) error->all(FLERR, "Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + + double k_one = utils::numeric(FLERR, arg[2], false, lmp); + double cut_one = utils::numeric(FLERR, arg[3], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { + k[i][j] = k_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairHarmonicCut::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); + cut[j][i] = cut[i][j]; + k[i][j] = mix_energy(k[i][i], k[j][j], cut[i][i], cut[j][j]); + k[j][i] = k[i][j]; + } + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i, j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j], sizeof(int), 1, fp); + if (setflag[i][j]) { + fwrite(&k[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairHarmonicCut::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i, j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); + } + MPI_Bcast(&k[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_restart_settings(FILE *fp) +{ + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairHarmonicCut::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); + } + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, k[i][i], cut[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairHarmonicCut::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) fprintf(fp, "%d %d %g %g\n", i, j, k[i][j], cut[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairHarmonicCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, double &fforce) +{ + if (rsq >= cutsq[itype][jtype]) { + fforce = 0.0; + return 0.0; + } + const double r = sqrt(rsq); + const double delta = cut[itype][jtype] - r; + const double philj = factor_lj * delta * k[itype][jtype]; + fforce = delta * philj / r; + return philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairHarmonicCut::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str, "k") == 0) return (void *) k; + if (strcmp(str, "cut") == 0) return (void *) cut; + return nullptr; +} diff --git a/src/EXTRA-PAIR/pair_harmonic_cut.h b/src/EXTRA-PAIR/pair_harmonic_cut.h new file mode 100644 index 0000000000..7dfdea995a --- /dev/null +++ b/src/EXTRA-PAIR/pair_harmonic_cut.h @@ -0,0 +1,72 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(harmonic/cut,PairHarmonicCut); +// clang-format on +#else + +#ifndef LMP_PAIR_HARMONIC_CUT_H +#define LMP_PAIR_HARMONIC_CUT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairHarmonicCut : public Pair { + public: + PairHarmonicCut(class LAMMPS *); + virtual ~PairHarmonicCut(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void write_data(FILE *); + void write_data_all(FILE *); + double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double **k, **cut; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +*/ diff --git a/src/OPENMP/pair_harmonic_cut_omp.cpp b/src/OPENMP/pair_harmonic_cut_omp.cpp new file mode 100644 index 0000000000..6fcdc50a63 --- /dev/null +++ b/src/OPENMP/pair_harmonic_cut_omp.cpp @@ -0,0 +1,152 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_harmonic_cut_omp.h" + +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neigh_list.h" +#include "suffix.h" + +#include "omp_compat.h" +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairHarmonicCutOMP::PairHarmonicCutOMP(LAMMPS *lmp) : PairHarmonicCut(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; +} + +/* ---------------------------------------------------------------------- */ + +void PairHarmonicCutOMP::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + if (evflag) { + if (eflag) { + if (force->newton_pair) + eval<1, 1, 1>(ifrom, ito, thr); + else + eval<1, 1, 0>(ifrom, ito, thr); + } else { + if (force->newton_pair) + eval<1, 0, 1>(ifrom, ito, thr); + else + eval<1, 0, 0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) + eval<0, 0, 1>(ifrom, ito, thr); + else + eval<0, 0, 0>(ifrom, ito, thr); + } + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +template +void PairHarmonicCutOMP::eval(int iifrom, int iito, ThrData *const thr) +{ + const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0]; + dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0]; + const int *_noalias const type = atom->type; + const double *_noalias const special_lj = force->special_lj; + const int *_noalias const ilist = list->ilist; + const int *_noalias const numneigh = list->numneigh; + const int *const *const firstneigh = list->firstneigh; + + double xtmp, ytmp, ztmp, delx, dely, delz, fxtmp, fytmp, fztmp; + double rsq, factor_lj; + + const int nlocal = atom->nlocal; + int j, jj, jnum, jtype; + + // loop over neighbors of my atoms + + for (int ii = iifrom; ii < iito; ++ii) { + const int i = ilist[ii]; + const int itype = type[i]; + const int *_noalias const jlist = firstneigh[i]; + const double *_noalias const cutsqi = cutsq[itype]; + + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + jnum = numneigh[i]; + fxtmp = fytmp = fztmp = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsqi[jtype]) { + const double r = sqrt(rsq); + const double delta = cut[itype][jtype] - r; + const double philj = factor_lj * delta * k[itype][jtype]; + const double fpair = delta * philj / r; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (NEWTON_PAIR || j < nlocal) { + f[j].x -= delx * fpair; + f[j].y -= dely * fpair; + f[j].z -= delz * fpair; + } + + if (EVFLAG) + ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, philj, 0.0, fpair, delx, dely, delz, thr); + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairHarmonicCutOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairHarmonicCut::memory_usage(); + + return bytes; +} diff --git a/src/OPENMP/pair_harmonic_cut_omp.h b/src/OPENMP/pair_harmonic_cut_omp.h new file mode 100644 index 0000000000..4f40d2cf2a --- /dev/null +++ b/src/OPENMP/pair_harmonic_cut_omp.h @@ -0,0 +1,48 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(harmonic/cut/omp,PairHarmonicCutOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_HARMONIC_CUT_OMP_H +#define LMP_PAIR_HARMONIC_CUT_OMP_H + +#include "pair_harmonic_cut.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairHarmonicCutOMP : public PairHarmonicCut, public ThrOMP { + + public: + PairHarmonicCutOMP(class LAMMPS *); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/unittest/force-styles/tests/mol-pair-harmonic_cut.yaml b/unittest/force-styles/tests/mol-pair-harmonic_cut.yaml new file mode 100644 index 0000000000..f1c39f1bd5 --- /dev/null +++ b/unittest/force-styles/tests/mol-pair-harmonic_cut.yaml @@ -0,0 +1,94 @@ +--- +lammps_version: 7 Jan 2022 +tags: generated +date_generated: Sat Jan 15 17:42:24 2022 +epsilon: 5e-14 +skip_tests: +prerequisites: ! | + atom full + pair harmonic/cut +pre_commands: ! "" +post_commands: ! | + pair_modify mix arithmetic +input_file: in.fourmol +pair_style: harmonic/cut +pair_coeff: ! | + 1 1 0.2 2.5 + 2 2 0.05 1.0 + 2 4 0.05 0.5 + 3 3 0.2 3.4 + 4 4 0.15 3.1 + 5 5 0.15 3.1 +extract: ! | + k 2 + cut 2 +natoms: 29 +init_vdwl: 0.6229303832094011 +init_coul: 0 +init_stress: ! |2- + 2.7086112336155066e-01 2.7116413086977192e-01 3.9219767031763964e-01 -1.0030847564854772e-01 4.4609473180383268e-02 -4.3637124113592815e-02 +init_forces: ! |2 + 1 8.0131849873788189e-03 4.9888010120265600e-02 2.8122360795900539e-02 + 2 7.9791227934915224e-03 6.5662021486410233e-03 -9.3912001601608869e-03 + 3 -3.0383030044476164e-02 -2.4462070024746502e-02 -6.7080526995988831e-03 + 4 -4.2871382862959282e-03 1.2944422156819623e-04 -3.1664690534881235e-03 + 5 -2.5653599407787658e-03 -5.6347847352330024e-03 5.7588483204582012e-03 + 6 -2.7432516621812987e-02 2.2443342050432843e-02 6.0842119146428050e-03 + 7 -2.3732010354850359e-02 4.8925005245008778e-03 -1.0554388509160603e-01 + 8 -1.9721449495438432e-02 8.7099489439665806e-03 7.7340656614040257e-02 + 9 3.1079843659201794e-03 3.3474898466584100e-03 1.3753907856580455e-02 + 10 1.5600383364335340e-02 -3.4154073795756010e-02 -2.4568763191891510e-02 + 11 -6.4786170866648325e-04 -2.7187839165446434e-03 -4.0027246972607657e-03 + 12 4.5446169193717613e-02 9.5334383719478614e-03 -2.3805265135058946e-02 + 13 3.6852167147282842e-03 -1.4578336695269954e-03 -6.9263908097426375e-05 + 14 -1.4516812720744084e-03 2.9539789066737132e-04 -3.7646771444256287e-03 + 15 -5.4837745225155548e-05 3.7892348244262858e-03 1.3279622443994443e-03 + 16 2.7804441490434852e-02 -2.9016569654227675e-02 -2.3977127174652516e-02 + 17 -1.3606174403879238e-03 -1.2150893147040220e-02 7.2609480510219010e-02 + 18 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 21 -2.8478436810585638e-03 -3.4997997565726423e-03 1.0311097670993996e-02 + 22 -5.0036838658335881e-03 -1.2125954307909393e-03 -7.8502223695393180e-03 + 23 7.8515275468921519e-03 4.7123951873635815e-03 -2.4608753014546784e-03 + 24 2.2422341594776187e-03 -9.2864513866737705e-03 5.4147040050256424e-03 + 25 -7.0931333559850684e-03 1.1418916865888746e-03 -5.9630640752421178e-03 + 26 4.8508991965074497e-03 8.1445597000848954e-03 5.4836007021647544e-04 + 27 1.4352981398486984e-03 -1.0223071626521522e-02 3.5886145449699289e-03 + 28 -7.4679157140535102e-03 3.2093693544577983e-03 -5.0528642356851916e-03 + 29 6.0326175742048118e-03 7.0137022720637250e-03 1.4642496907152625e-03 +run_vdwl: 0.6227892941470253 +run_coul: 0 +run_stress: ! |2- + 2.7079595899259656e-01 2.7106155255810421e-01 3.9199181631609364e-01 -1.0017597886837336e-01 4.4639512685950337e-02 -4.3538044335882868e-02 +run_forces: ! |2 + 1 8.0241133736762959e-03 4.9859217970288555e-02 2.8099645116503875e-02 + 2 7.9858506322803297e-03 6.5851530052680168e-03 -9.3657942914726639e-03 + 3 -3.0378796721416249e-02 -2.4495674709975927e-02 -6.7374520329115941e-03 + 4 -4.2659306402640165e-03 1.3406094892709258e-04 -3.1564081120146015e-03 + 5 -2.5644106889098472e-03 -5.6291258903991056e-03 5.7537772274397249e-03 + 6 -2.7398860728585855e-02 2.2446081967441193e-02 6.0769742046490674e-03 + 7 -2.3725250988273959e-02 4.8748360958146379e-03 -1.0543946577300001e-01 + 8 -1.9795760001146530e-02 8.7619549352892182e-03 7.7257651402956734e-02 + 9 3.1107511049119794e-03 3.3343639112369677e-03 1.3749677634153570e-02 + 10 1.5635110483272714e-02 -3.4177809303542570e-02 -2.4562568959156386e-02 + 11 -6.4680122415666373e-04 -2.7045376093154955e-03 -3.9833439702283155e-03 + 12 4.5410386808201031e-02 9.5232810722562407e-03 -2.3816427499003763e-02 + 13 3.6819606745360318e-03 -1.4490090012638590e-03 -6.8397423783681964e-05 + 14 -1.4472530230181654e-03 2.8919520744570527e-04 -3.7524156545238243e-03 + 15 -5.6155855006963267e-05 3.7938646788886268e-03 1.3339639717048200e-03 + 16 2.7787095304862428e-02 -2.9020839215424241e-02 -2.3998357673493245e-02 + 17 -1.3560485109625556e-03 -1.2125014062935051e-02 7.2608941832180307e-02 + 18 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 21 -2.8538862991433229e-03 -3.4800313716197438e-03 1.0305823645864499e-02 + 22 -5.0214625596301642e-03 -1.2260879074991976e-03 -7.8500736428379052e-03 + 23 7.8753488587734871e-03 4.7061192791189414e-03 -2.4557500030265946e-03 + 24 2.2600769444580133e-03 -9.2990099176050046e-03 5.4253104994500942e-03 + 25 -7.1247130587123505e-03 1.1308496354809658e-03 -5.9913273289481242e-03 + 26 4.8646361142543372e-03 8.1681602821240384e-03 5.6601682949803040e-04 + 27 1.4504263375229804e-03 -1.0234012017965840e-02 3.5742943025157811e-03 + 28 -7.4838995855112476e-03 3.2116101208619755e-03 -5.0534250649872060e-03 + 29 6.0334732479882672e-03 7.0224018971038644e-03 1.4791307624714251e-03 +... From 943fe487b54f165796c77a2e1b2a6b44012f3a9d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 16 Jan 2022 15:36:01 -0500 Subject: [PATCH 11/75] update whitespace and argument formats for longer source lines --- src/molecule.cpp | 222 +++++++++++++++++------------------------------ 1 file changed, 81 insertions(+), 141 deletions(-) diff --git a/src/molecule.cpp b/src/molecule.cpp index 40ca218ecc..d6c839dfc4 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -57,8 +57,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : id = utils::strdup(arg[0]); if (!utils::is_id(id)) - error->all(FLERR,"Molecule template ID must have only " - "alphanumeric or underscore characters"); + error->all(FLERR,"Molecule template ID must have only alphanumeric or underscore characters"); // parse args until reach unknown arg (next file) @@ -130,8 +129,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : if (me == 0) { fp = fopen(arg[ifile],"r"); if (fp == nullptr) - error->one(FLERR,"Cannot open molecule file {}: {}", - arg[ifile], utils::getsyserror()); + error->one(FLERR,"Cannot open molecule file {}: {}",arg[ifile], utils::getsyserror()); } read(0); if (me == 0) fclose(fp); @@ -496,8 +494,7 @@ void Molecule::read(int flag) if (nmatch != nwant) error->one(FLERR,"Invalid header line format in molecule file"); } catch (TokenizerException &e) { - error->one(FLERR, "Invalid header in molecule file\n" - "{}", e.what()); + error->one(FLERR,"Invalid header in molecule file: {}", e.what()); } } @@ -618,11 +615,9 @@ void Molecule::read(int flag) // Error: Either a too long/short section or a typo in the keyword if (utils::strmatch(keyword,"^[A-Za-z ]+$")) - error->one(FLERR,"Unknown section '{}' in molecule " - "file\n",keyword); + error->one(FLERR,"Unknown section '{}' in molecule file\n",keyword); else error->one(FLERR,"Unexpected line in molecule file " - "while looking for the next " - "section:\n{}",line); + "while looking for the next section:\n{}",line); } keyword = parse_keyword(1,line); } @@ -648,8 +643,7 @@ void Molecule::read(int flag) if (bondflag && specialflag == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot auto-generate special bonds before " - "simulation box is defined"); + error->all(FLERR,"Cannot auto-generate special bonds before simulation box is defined"); if (flag) { special_generate(); @@ -691,8 +685,7 @@ void Molecule::coords(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 4) - error->all(FLERR,"Invalid line in Coords section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Coords section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -707,19 +700,16 @@ void Molecule::coords(char *line) x[iatom][2] *= sizescale; } } catch (TokenizerException &e) { - error->all(FLERR,"Invalid line in Coords section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Coords section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Coords " - "section of molecule file",i+1); + if (count[i] == 0) error->all(FLERR,"Atom {} missing in Coords section of molecule file",i+1); if (domain->dimension == 2) { for (int i = 0; i < natoms; i++) if (x[i][2] != 0.0) - error->all(FLERR,"Z coord in molecule file for atom {} " - "must be 0.0 for 2d-simulation.",i+1); + error->all(FLERR,"Z coord in molecule file for atom {} must be 0.0 for 2d-simulation",i+1); } } @@ -737,8 +727,7 @@ void Molecule::types(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Types section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Types section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -748,17 +737,14 @@ void Molecule::types(char *line) type[iatom] += toffset; } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Types section of " - "molecule file: {}\n{}", e.what(),line); + error->all(FLERR,"Invalid line in Types section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) { - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Types " - "section of molecule file",i+1); + if (count[i] == 0) error->all(FLERR,"Atom {} missing in Types section of molecule file",i+1); if ((type[i] <= 0) || (domain->box_exist && (type[i] > atom->ntypes))) - error->all(FLERR,"Invalid atom type {} for atom {} " - "in molecule file",type[i],i+1); + error->all(FLERR,"Invalid atom type {} for atom {} in molecule file",type[i],i+1); ntypes = MAX(ntypes,type[i]); } @@ -777,8 +763,7 @@ void Molecule::molecules(char *line) readline(line); ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Molecules section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Molecules section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -788,19 +773,17 @@ void Molecule::molecules(char *line) // molecule[iatom] += moffset; // placeholder for possible molecule offset } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Molecules section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Molecules section of molecule file: {}\n{}",e.what(),line); } - for (int i = 0; i < natoms; i++) - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Molecules " - "section of molecule file",i+1); - - for (int i = 0; i < natoms; i++) + for (int i = 0; i < natoms; i++) { + if (count[i] == 0) + error->all(FLERR,"Atom {} missing in Molecules section of molecule file",i+1); + } + for (int i = 0; i < natoms; i++) { if (molecule[i] < 0) - error->all(FLERR,"Invalid molecule ID {} for atom {} " - "in molecule file",molecule[i],i+1); - + error->all(FLERR,"Invalid molecule ID {} for atom {} in molecule file",molecule[i],i+1); + } for (int i = 0; i < natoms; i++) nmolecules = MAX(nmolecules,molecule[i]); } @@ -818,23 +801,21 @@ void Molecule::fragments(char *line) ValueTokenizer values(utils::trim_comment(line)); if ((int)values.count() > natoms+1) - error->all(FLERR,"Too many atoms per fragment in Fragments " - "section of molecule file"); + error->all(FLERR,"Too many atoms per fragment in Fragments section of molecule file"); fragmentnames[i] = values.next_string(); while (values.has_next()) { int iatom = values.next_int()-1; if (iatom < 0 || iatom >= natoms) - error->all(FLERR,"Invalid atom ID {} for fragment {} in " - "Fragments section of molecule file", - iatom+1, fragmentnames[i]); + error->all(FLERR,"Invalid atom ID {} for fragment {} in Fragments section of " + "molecule file", iatom+1, fragmentnames[i]); fragmentmask[i][iatom] = 1; } } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid atom ID in Fragments section of " - "molecule file: {}\n{}", e.what(),line); + error->all(FLERR,"Invalid atom ID in Fragments section of " + "molecule file: {}\n{}", e.what(),line); } } @@ -851,8 +832,7 @@ void Molecule::charges(char *line) ValueTokenizer values(utils::trim_comment(line)); if ((int)values.count() != 2) - error->all(FLERR,"Invalid line in Charges section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Charges section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -862,13 +842,13 @@ void Molecule::charges(char *line) q[iatom] = values.next_double(); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Charges section of " - "molecule file: {}.\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Charges section of molecule file: {}\n{}",e.what(),line); } - for (int i = 0; i < natoms; i++) - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Charges " - "section of molecule file",i+1); + for (int i = 0; i < natoms; i++) { + if (count[i] == 0) + error->all(FLERR,"Atom {} missing in Charges section of molecule file",i+1); + } } /* ---------------------------------------------------------------------- @@ -885,8 +865,7 @@ void Molecule::diameters(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Diameters section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Diameters section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) error->all(FLERR,"Invalid atom index in Diameters section of molecule file"); @@ -897,16 +876,14 @@ void Molecule::diameters(char *line) maxradius = MAX(maxradius,radius[iatom]); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Diameters section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Diameters section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) { - if (count[i] == 0) error->all(FLERR,"Atom {} missing in Diameters " - "section of molecule file",i+1); + if (count[i] == 0) + error->all(FLERR,"Atom {} missing in Diameters section of molecule file",i+1); if (radius[i] < 0.0) - error->all(FLERR,"Invalid atom diameter {} for atom {} " - "in molecule file", radius[i], i+1); + error->all(FLERR,"Invalid atom diameter {} for atom {} in molecule file", radius[i], i+1); } } @@ -923,8 +900,7 @@ void Molecule::masses(char *line) ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 2) - error->all(FLERR,"Invalid line in Masses section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Masses section of molecule file: {}",line); int iatom = values.next_int() - 1; if (iatom < 0 || iatom >= natoms) @@ -934,8 +910,7 @@ void Molecule::masses(char *line) rmass[iatom] *= sizescale*sizescale*sizescale; } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Masses section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Masses section of molecule file: {}\n{}",e.what(),line); } for (int i = 0; i < natoms; i++) { @@ -972,15 +947,13 @@ void Molecule::bonds(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 4) - error->all(FLERR,"Invalid line in Bonds section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Bonds section of molecule file: {}",line); values.next_int(); itype = values.next_int(); atom1 = values.next_tagint(); atom2 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Bonds section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Bonds section of molecule file: {}\n{}",e.what(),line); } itype += boffset; @@ -1042,16 +1015,14 @@ void Molecule::angles(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 5) - error->all(FLERR,"Invalid line in Angles section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Angles section of molecule file: {}",line); values.next_int(); itype = values.next_int(); atom1 = values.next_tagint(); atom2 = values.next_tagint(); atom3 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Angles section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Angles section of molecule file: {}\n{}",e.what(),line); } itype += aoffset; @@ -1128,8 +1099,7 @@ void Molecule::dihedrals(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 6) - error->all(FLERR,"Invalid line in Dihedrals section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Dihedrals section of molecule file: {}",line); values.next_int(); itype = values.next_int(); @@ -1138,8 +1108,7 @@ void Molecule::dihedrals(int flag, char *line) atom3 = values.next_tagint(); atom4 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Dihedrals section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Dihedrals section of molecule file: {}\n{}",e.what(),line); } itype += doffset; @@ -1150,8 +1119,7 @@ void Molecule::dihedrals(int flag, char *line) (atom4 <= 0) || (atom4 > natoms) || (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) - error->all(FLERR, - "Invalid atom ID in dihedrals section of molecule file"); + error->all(FLERR,"Invalid atom ID in dihedrals section of molecule file"); if ((itype <= 0) || (domain->box_exist && (itype > atom->ndihedraltypes))) error->all(FLERR,"Invalid dihedral type in Dihedrals section of molecule file"); @@ -1230,8 +1198,7 @@ void Molecule::impropers(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 6) - error->all(FLERR,"Invalid line in Impropers section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Impropers section of molecule file: {}",line); values.next_int(); itype = values.next_int(); atom1 = values.next_tagint(); @@ -1239,8 +1206,7 @@ void Molecule::impropers(int flag, char *line) atom3 = values.next_tagint(); atom4 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Impropers section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Impropers section of molecule file: {}\n{}",e.what(),line); } itype += ioffset; @@ -1251,8 +1217,7 @@ void Molecule::impropers(int flag, char *line) (atom4 <= 0) || (atom4 > natoms) || (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) - error->all(FLERR, - "Invalid atom ID in impropers section of molecule file"); + error->all(FLERR,"Invalid atom ID in impropers section of molecule file"); if ((itype <= 0) || (domain->box_exist && (itype > atom->nimpropertypes))) error->all(FLERR,"Invalid improper type in Impropers section of molecule file"); @@ -1325,15 +1290,14 @@ void Molecule::nspecial_read(int flag, char *line) try { ValueTokenizer values(utils::trim_comment(line)); if (values.count() != 4) - error->all(FLERR,"Invalid line in Special Bond Counts section of " - "molecule file: {}",line); + error->all(FLERR,"Invalid line in Special Bond Counts section of molecule file: {}",line); values.next_int(); c1 = values.next_tagint(); c2 = values.next_tagint(); c3 = values.next_tagint(); } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Special Bond Counts section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Special Bond Counts section of " + "molecule file: {}\n{}",e.what(),line); } if (flag) { @@ -1358,8 +1322,7 @@ void Molecule::special_read(char *line) int nwords = values.count(); if (nwords != nspecial[i][2]+1) - error->all(FLERR,"Molecule file special list " - "does not match special count"); + error->all(FLERR,"Molecule file special list does not match special count"); values.next_int(); // ignore @@ -1367,13 +1330,12 @@ void Molecule::special_read(char *line) special[i][m-1] = values.next_tagint(); if (special[i][m-1] <= 0 || special[i][m-1] > natoms || special[i][m-1] == i+1) - error->all(FLERR,"Invalid atom index in Special Bonds " - "section of molecule file"); + error->all(FLERR,"Invalid atom index in Special Bonds section of molecule file"); } } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid line in Special Bonds section of " - "molecule file: {}\n{}",e.what(),line); + error->all(FLERR,"Invalid line in Special Bonds section of " + "molecule file: {}\n{}",e.what(),line); } } @@ -1502,8 +1464,7 @@ void Molecule::shakeflag_read(char *line) shake_flag[i] = values.next_int(); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid Shake Flags section in molecule file\n" - "{}", e.what()); + error->all(FLERR,"Invalid Shake Flags section in molecule file: {}", e.what()); } for (int i = 0; i < natoms; i++) @@ -1572,8 +1533,7 @@ void Molecule::shakeatom_read(char *line) } } catch (TokenizerException &e) { - error->all(FLERR,"Invalid shake atom in molecule file\n" - "{}", e.what()); + error->all(FLERR,"Invalid shake atom in molecule file: {}", e.what()); } for (int i = 0; i < natoms; i++) { @@ -1642,8 +1602,7 @@ void Molecule::shaketype_read(char *line) error->all(FLERR,"Invalid shake type data in molecule file"); } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid shake type data in molecule file\n", - "{}", e.what()); + error->all(FLERR,"Invalid shake type data in molecule file: {}",e.what()); } for (int i = 0; i < natoms; i++) { @@ -1695,8 +1654,7 @@ void Molecule::body(int flag, int pflag, char *line) } else nword += ncount; } } catch (TokenizerException &e) { - error->all(FLERR, "Invalid body params in molecule file\n", - "{}", e.what()); + error->all(FLERR,"Invalid body params in molecule file: {}", e.what()); } } @@ -1735,8 +1693,7 @@ void Molecule::check_attributes(int flag) if (onemol->rmassflag && !atom->rmass_flag) mismatch = 1; if (mismatch && me == 0) - error->warning(FLERR, - "Molecule attributes do not match system attributes"); + error->warning(FLERR,"Molecule attributes do not match system attributes"); // for all atom styles, check nbondtype,etc @@ -1770,8 +1727,7 @@ void Molecule::check_attributes(int flag) // warn if molecule topology defined but no special settings if (onemol->bondflag && !onemol->specialflag) - if (me == 0) error->warning(FLERR,"Molecule has bond topology " - "but no special bond settings"); + if (me == 0) error->warning(FLERR,"Molecule has bond topology but no special bond settings"); } } @@ -1878,47 +1834,31 @@ void Molecule::allocate() memory->create(special,natoms,maxspecial,"molecule:special"); if (bondflag) { - memory->create(bond_type,natoms,bond_per_atom, - "molecule:bond_type"); - memory->create(bond_atom,natoms,bond_per_atom, - "molecule:bond_atom"); + memory->create(bond_type,natoms,bond_per_atom,"molecule:bond_type"); + memory->create(bond_atom,natoms,bond_per_atom,"molecule:bond_atom"); } if (angleflag) { - memory->create(angle_type,natoms,angle_per_atom, - "molecule:angle_type"); - memory->create(angle_atom1,natoms,angle_per_atom, - "molecule:angle_atom1"); - memory->create(angle_atom2,natoms,angle_per_atom, - "molecule:angle_atom2"); - memory->create(angle_atom3,natoms,angle_per_atom, - "molecule:angle_atom3"); + memory->create(angle_type,natoms,angle_per_atom,"molecule:angle_type"); + memory->create(angle_atom1,natoms,angle_per_atom,"molecule:angle_atom1"); + memory->create(angle_atom2,natoms,angle_per_atom,"molecule:angle_atom2"); + memory->create(angle_atom3,natoms,angle_per_atom,"molecule:angle_atom3"); } if (dihedralflag) { - memory->create(dihedral_type,natoms,dihedral_per_atom, - "molecule:dihedral_type"); - memory->create(dihedral_atom1,natoms,dihedral_per_atom, - "molecule:dihedral_atom1"); - memory->create(dihedral_atom2,natoms,dihedral_per_atom, - "molecule:dihedral_atom2"); - memory->create(dihedral_atom3,natoms,dihedral_per_atom, - "molecule:dihedral_atom3"); - memory->create(dihedral_atom4,natoms,dihedral_per_atom, - "molecule:dihedral_atom4"); + memory->create(dihedral_type,natoms,dihedral_per_atom,"molecule:dihedral_type"); + memory->create(dihedral_atom1,natoms,dihedral_per_atom,"molecule:dihedral_atom1"); + memory->create(dihedral_atom2,natoms,dihedral_per_atom,"molecule:dihedral_atom2"); + memory->create(dihedral_atom3,natoms,dihedral_per_atom,"molecule:dihedral_atom3"); + memory->create(dihedral_atom4,natoms,dihedral_per_atom,"molecule:dihedral_atom4"); } if (improperflag) { - memory->create(improper_type,natoms,improper_per_atom, - "molecule:improper_type"); - memory->create(improper_atom1,natoms,improper_per_atom, - "molecule:improper_atom1"); - memory->create(improper_atom2,natoms,improper_per_atom, - "molecule:improper_atom2"); - memory->create(improper_atom3,natoms,improper_per_atom, - "molecule:improper_atom3"); - memory->create(improper_atom4,natoms,improper_per_atom, - "molecule:improper_atom4"); + memory->create(improper_type,natoms,improper_per_atom,"molecule:improper_type"); + memory->create(improper_atom1,natoms,improper_per_atom,"molecule:improper_atom1"); + memory->create(improper_atom2,natoms,improper_per_atom,"molecule:improper_atom2"); + memory->create(improper_atom3,natoms,improper_per_atom,"molecule:improper_atom3"); + memory->create(improper_atom4,natoms,improper_per_atom,"molecule:improper_atom4"); } if (shakeflag) { @@ -2058,7 +1998,7 @@ void Molecule::skip_lines(int n, char *line, const std::string §ion) readline(line); if (utils::strmatch(utils::trim(utils::trim_comment(line)),"^[A-Za-z ]+$")) error->one(FLERR,"Unexpected line in molecule file while " - "skipping {} section:\n{}",section,line); + "skipping {} section:\n{}",section,line); } } From 0eeb3b203cf84fcde15a6b34f2c239d8045ba33d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 16 Jan 2022 16:50:23 -0500 Subject: [PATCH 12/75] add tests for molecule command --- unittest/formats/test_molecule_file.cpp | 27 +++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/unittest/formats/test_molecule_file.cpp b/unittest/formats/test_molecule_file.cpp index 2cca7a6832..e802ebfa2c 100644 --- a/unittest/formats/test_molecule_file.cpp +++ b/unittest/formats/test_molecule_file.cpp @@ -176,6 +176,33 @@ TEST_F(MoleculeFileTest, minimal) ASSERT_THAT(output, MatchesRegex(".*Read molecule template.*1 molecules.*1 atoms.*0 bonds.*")); } +TEST_F(MoleculeFileTest, notype) +{ + BEGIN_CAPTURE_OUTPUT(); + command("atom_style atomic"); + command("region box block 0 1 0 1 0 1"); + command("create_box 1 box"); + run_mol_cmd(test_name, "", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"); + auto output = END_CAPTURE_OUTPUT(); + ASSERT_THAT(output, MatchesRegex(".*Read molecule template.*1 molecules.*1 atoms.*0 bonds.*")); + TEST_FAILURE(".*ERROR: Create_atoms molecule must have atom types.*", + command("create_atoms 0 single 0.0 0.0 0.0 mol notype 542465");); + } + +TEST_F(MoleculeFileTest, extramass) +{ + BEGIN_CAPTURE_OUTPUT(); + command("atom_style atomic"); + command("region box block 0 1 0 1 0 1"); + command("create_box 1 box"); + run_mol_cmd(test_name, "", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n" + " Types\n\n 1 1\n Masses\n\n 1 1.0\n"); + command("create_atoms 0 single 0.0 0.0 0.0 mol extramass 73546"); + auto output = END_CAPTURE_OUTPUT(); + ASSERT_THAT(output, MatchesRegex(".*WARNING: Molecule attributes do not match " + "system attributes.*")); +} + TEST_F(MoleculeFileTest, twomols) { BEGIN_CAPTURE_OUTPUT(); From dc6e558191c756b64e9e01fc7505a40328c54491 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 16 Jan 2022 20:20:07 -0500 Subject: [PATCH 13/75] use Tokenizer class to parse bond colors --- src/dump_image.cpp | 27 +++++++-------------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/src/dump_image.cpp b/src/dump_image.cpp index c073d152f8..9aecb58120 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -1200,11 +1200,11 @@ int DumpImage::modify_param(int narg, char **arg) utils::bounds(FLERR,arg[1],1,atom->ntypes,nlo,nhi,error); // get list of colors + // assign colors in round-robin fashion to types + auto colors = Tokenizer(arg[2],"/").as_vector(); const int ncolors = colors.size(); - // assign colors in round-robin fashion to types - int m = 0; for (int i = nlo; i <= nhi; i++) { colortype[i] = image->color2rgb(colors[m%ncolors].c_str()); @@ -1249,32 +1249,19 @@ int DumpImage::modify_param(int narg, char **arg) int nlo,nhi; utils::bounds(FLERR,arg[1],1,atom->nbondtypes,nlo,nhi,error); - // ptrs = list of ncount colornames separated by '/' + // process list of ncount colornames separated by '/' + // assign colors in round-robin fashion to bond types - int ncount = 1; - char *nextptr; - char *ptr = arg[2]; - while ((nextptr = strchr(ptr,'/'))) { - ptr = nextptr + 1; - ncount++; - } - char **ptrs = new char*[ncount+1]; - ncount = 0; - ptrs[ncount++] = strtok(arg[2],"/"); - while ((ptrs[ncount++] = strtok(nullptr,"/"))); - ncount--; - - // assign each of ncount colors in round-robin fashion to types + auto colors = Tokenizer(arg[2],"/").as_vector(); + const int ncolors = colors.size(); int m = 0; for (int i = nlo; i <= nhi; i++) { - bcolortype[i] = image->color2rgb(ptrs[m%ncount]); + bcolortype[i] = image->color2rgb(colors[m%ncolors].c_str()); if (bcolortype[i] == nullptr) error->all(FLERR,"Invalid color in dump_modify command"); m++; } - - delete [] ptrs; return 3; } From cb796e8b601f0a29918a864fa84bf9daf8134451 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 17 Jan 2022 08:21:52 -0700 Subject: [PATCH 14/75] Fix HIP compile issues --- lib/kokkos/Makefile.kokkos | 3 +++ lib/kokkos/Makefile.targets | 2 ++ src/KOKKOS/fix_acks2_reaxff_kokkos.cpp | 2 ++ 3 files changed, 7 insertions(+) diff --git a/lib/kokkos/Makefile.kokkos b/lib/kokkos/Makefile.kokkos index 7ffea5a62c..899d6e49a2 100644 --- a/lib/kokkos/Makefile.kokkos +++ b/lib/kokkos/Makefile.kokkos @@ -1224,6 +1224,9 @@ ifeq ($(KOKKOS_INTERNAL_USE_HIP), 1) KOKKOS_SRC += $(wildcard $(KOKKOS_PATH)/core/src/HIP/*.cpp) + ifeq ($(KOKKOS_INTERNAL_ENABLE_DESUL_ATOMICS), 1) + KOKKOS_SRC += $(KOKKOS_PATH)/core/src/desul/src/Lock_Array_HIP.cpp + endif KOKKOS_HEADERS += $(wildcard $(KOKKOS_PATH)/core/src/HIP/*.hpp) KOKKOS_CXXFLAGS+=$(KOKKOS_INTERNAL_HIP_ARCH_FLAG) diff --git a/lib/kokkos/Makefile.targets b/lib/kokkos/Makefile.targets index 93854d0cf1..c097e80fec 100644 --- a/lib/kokkos/Makefile.targets +++ b/lib/kokkos/Makefile.targets @@ -68,6 +68,8 @@ Kokkos_HIP_Instance.o: $(KOKKOS_CPP_DEPENDS) $(KOKKOS_PATH)/core/src/HIP/Kokkos_ $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) -c $(KOKKOS_PATH)/core/src/HIP/Kokkos_HIP_Instance.cpp Kokkos_HIP_Locks.o: $(KOKKOS_CPP_DEPENDS) $(KOKKOS_PATH)/core/src/HIP/Kokkos_HIP_Locks.cpp $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) -c $(KOKKOS_PATH)/core/src/HIP/Kokkos_HIP_Locks.cpp +Lock_Array_HIP.o: $(KOKKOS_CPP_DEPENDS) $(KOKKOS_PATH)/core/src/desul/src/Lock_Array_HIP.cpp + $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) -c $(KOKKOS_PATH)/core/src/desul/src/Lock_Array_HIP.cpp endif ifeq ($(KOKKOS_INTERNAL_USE_PTHREADS), 1) diff --git a/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp b/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp index 8379dc8f46..582e9a39ce 100644 --- a/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp +++ b/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp @@ -633,6 +633,7 @@ void FixACKS2ReaxFFKokkos::compute_h_item(int ii, int &m_fill, const template template +KOKKOS_INLINE_FUNCTION void FixACKS2ReaxFFKokkos::compute_h_team( const typename Kokkos::TeamPolicy::member_type &team, int atoms_per_team, int vector_length) const { @@ -935,6 +936,7 @@ void FixACKS2ReaxFFKokkos::compute_x_item(int ii, int &m_fill, const template template +KOKKOS_INLINE_FUNCTION void FixACKS2ReaxFFKokkos::compute_x_team( const typename Kokkos::TeamPolicy::member_type &team, int atoms_per_team, int vector_length) const { From a93e5baa73d3f108f2503bd7cfcd41d252c1b115 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 17 Jan 2022 08:44:49 -0700 Subject: [PATCH 15/75] Add Kokkos HIP Makefiles --- src/MAKE/MACHINES/Makefile.crusher_kokkos | 123 ++++++++++++++++++++++ src/MAKE/MACHINES/Makefile.spock_kokkos | 123 ++++++++++++++++++++++ 2 files changed, 246 insertions(+) create mode 100644 src/MAKE/MACHINES/Makefile.crusher_kokkos create mode 100644 src/MAKE/MACHINES/Makefile.spock_kokkos diff --git a/src/MAKE/MACHINES/Makefile.crusher_kokkos b/src/MAKE/MACHINES/Makefile.crusher_kokkos new file mode 100644 index 0000000000..167341ec9c --- /dev/null +++ b/src/MAKE/MACHINES/Makefile.crusher_kokkos @@ -0,0 +1,123 @@ +# summit_kokkos = KOKKOS/CUDA, V100 GPU and Power9, IBM Spectrum MPI, nvcc compiler with gcc + +SHELL = /bin/sh + +# --------------------------------------------------------------------- +# compiler/linker settings +# specify flags and libraries needed for your compiler + +KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) + +CC = hipcc +CCFLAGS = -g -O3 -munsafe-fp-atomics -DNDEBUG +SHFLAGS = -fPIC +DEPFLAGS = -M + +LINK = hipcc +LINKFLAGS = -g -O3 +LIB = +SIZE = size + +ARCHIVE = ar +ARFLAGS = -rc +SHLIBFLAGS = -shared +KOKKOS_DEVICES = HIP +KOKKOS_ARCH = Zen3,Vega90A + +# --------------------------------------------------------------------- +# LAMMPS-specific settings, all OPTIONAL +# specify settings for LAMMPS features you will use +# if you change any -D setting, do full re-compile after "make clean" + +# LAMMPS ifdef settings +# see possible settings in Section 3.5 of the manual + +LMP_INC = -DLAMMPS_GZIP + +# MPI library +# see discussion in Section 3.4 of the manual +# MPI wrapper compiler/linker can provide this info +# can point to dummy MPI library in src/STUBS as in Makefile.serial +# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts +# INC = path for mpi.h, MPI compiler settings +# PATH = path for MPI library +# LIB = name of MPI library + +MPI_INC = -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1 -I${MPICH_DIR}/include +MPI_PATH = +MPI_LIB = -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa + +# FFT library +# see discussion in Section 3.5.2 of manual +# can be left blank to use provided KISS FFT library +# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings +# PATH = path for FFT library +# LIB = name of FFT library + +FFT_INC = +FFT_PATH = +FFT_LIB = + +# JPEG and/or PNG library +# see discussion in Section 3.5.4 of manual +# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC +# INC = path(s) for jpeglib.h and/or png.h +# PATH = path(s) for JPEG library and/or PNG library +# LIB = name(s) of JPEG library and/or PNG library + +JPG_INC = +JPG_PATH = +JPG_LIB = + +# --------------------------------------------------------------------- +# build rules and dependencies +# do not edit this section + +include Makefile.package.settings +include Makefile.package + +EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC) +EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH) +EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB) +EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS) +EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS) + +# Path to src files + +vpath %.cpp .. +vpath %.h .. + +# Link target + +$(EXE): main.o $(LMPLIB) $(EXTRA_LINK_DEPENDS) + $(LINK) $(LINKFLAGS) main.o $(EXTRA_PATH) $(LMPLINK) $(EXTRA_LIB) $(LIB) -o $@ + $(SIZE) $@ + +# Library targets + +$(ARLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) + @rm -f ../$(ARLIB) + $(ARCHIVE) $(ARFLAGS) ../$(ARLIB) $(OBJ) + @rm -f $(ARLIB) + @ln -s ../$(ARLIB) $(ARLIB) + +$(SHLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) + $(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o ../$(SHLIB) \ + $(OBJ) $(EXTRA_LIB) $(LIB) + @rm -f $(SHLIB) + @ln -s ../$(SHLIB) $(SHLIB) + +# Compilation rules + +%.o:%.cpp + $(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $< + +# Individual dependencies + +depend : fastdep.exe $(SRC) + @./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1 + +fastdep.exe: ../DEPEND/fastdep.c + cc -O -o $@ $< + +sinclude .depend diff --git a/src/MAKE/MACHINES/Makefile.spock_kokkos b/src/MAKE/MACHINES/Makefile.spock_kokkos new file mode 100644 index 0000000000..44fd213962 --- /dev/null +++ b/src/MAKE/MACHINES/Makefile.spock_kokkos @@ -0,0 +1,123 @@ +# summit_kokkos = KOKKOS/CUDA, V100 GPU and Power9, IBM Spectrum MPI, nvcc compiler with gcc + +SHELL = /bin/sh + +# --------------------------------------------------------------------- +# compiler/linker settings +# specify flags and libraries needed for your compiler + +KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) + +CC = hipcc +CCFLAGS = -g -O3 -DNDEBUG +SHFLAGS = -fPIC +DEPFLAGS = -M + +LINK = hipcc +LINKFLAGS = -g -O3 +LIB = +SIZE = size + +ARCHIVE = ar +ARFLAGS = -rc +SHLIBFLAGS = -shared +KOKKOS_DEVICES = HIP +KOKKOS_ARCH = Zen2,Vega908 + +# --------------------------------------------------------------------- +# LAMMPS-specific settings, all OPTIONAL +# specify settings for LAMMPS features you will use +# if you change any -D setting, do full re-compile after "make clean" + +# LAMMPS ifdef settings +# see possible settings in Section 3.5 of the manual + +LMP_INC = -DLAMMPS_GZIP + +# MPI library +# see discussion in Section 3.4 of the manual +# MPI wrapper compiler/linker can provide this info +# can point to dummy MPI library in src/STUBS as in Makefile.serial +# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts +# INC = path for mpi.h, MPI compiler settings +# PATH = path for MPI library +# LIB = name of MPI library + +MPI_INC = -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1 -I${MPICH_DIR}/include +MPI_PATH = +MPI_LIB = -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa + +# FFT library +# see discussion in Section 3.5.2 of manual +# can be left blank to use provided KISS FFT library +# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings +# PATH = path for FFT library +# LIB = name of FFT library + +FFT_INC = +FFT_PATH = +FFT_LIB = + +# JPEG and/or PNG library +# see discussion in Section 3.5.4 of manual +# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC +# INC = path(s) for jpeglib.h and/or png.h +# PATH = path(s) for JPEG library and/or PNG library +# LIB = name(s) of JPEG library and/or PNG library + +JPG_INC = +JPG_PATH = +JPG_LIB = + +# --------------------------------------------------------------------- +# build rules and dependencies +# do not edit this section + +include Makefile.package.settings +include Makefile.package + +EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC) +EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH) +EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB) +EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS) +EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS) + +# Path to src files + +vpath %.cpp .. +vpath %.h .. + +# Link target + +$(EXE): main.o $(LMPLIB) $(EXTRA_LINK_DEPENDS) + $(LINK) $(LINKFLAGS) main.o $(EXTRA_PATH) $(LMPLINK) $(EXTRA_LIB) $(LIB) -o $@ + $(SIZE) $@ + +# Library targets + +$(ARLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) + @rm -f ../$(ARLIB) + $(ARCHIVE) $(ARFLAGS) ../$(ARLIB) $(OBJ) + @rm -f $(ARLIB) + @ln -s ../$(ARLIB) $(ARLIB) + +$(SHLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) + $(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o ../$(SHLIB) \ + $(OBJ) $(EXTRA_LIB) $(LIB) + @rm -f $(SHLIB) + @ln -s ../$(SHLIB) $(SHLIB) + +# Compilation rules + +%.o:%.cpp + $(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $< + +# Individual dependencies + +depend : fastdep.exe $(SRC) + @./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1 + +fastdep.exe: ../DEPEND/fastdep.c + cc -O -o $@ $< + +sinclude .depend From 8b89be60619ba8879fee915ce4d01a93f4574d39 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 17 Jan 2022 08:49:00 -0700 Subject: [PATCH 16/75] Kokkos SNAP tuning for HIP --- src/KOKKOS/pair_snap_kokkos.h | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/KOKKOS/pair_snap_kokkos.h b/src/KOKKOS/pair_snap_kokkos.h index bd56d87f59..e7af536782 100644 --- a/src/KOKKOS/pair_snap_kokkos.h +++ b/src/KOKKOS/pair_snap_kokkos.h @@ -85,6 +85,19 @@ public: using complex = SNAComplex; // Static team/tile sizes for device offload + +#ifdef KOKKOS_ENABLE_HIP + static constexpr int team_size_compute_neigh = 2; + static constexpr int tile_size_compute_ck = 2; + static constexpr int tile_size_pre_ui = 2; + static constexpr int team_size_compute_ui = 2; + static constexpr int tile_size_transform_ui = 2; + static constexpr int tile_size_compute_zi = 2; + static constexpr int tile_size_compute_bi = 2; + static constexpr int tile_size_transform_bi = 2; + static constexpr int tile_size_compute_yi = 2; + static constexpr int team_size_compute_fused_deidrj = 2; +#else static constexpr int team_size_compute_neigh = 4; static constexpr int tile_size_compute_ck = 4; static constexpr int tile_size_pre_ui = 4; @@ -95,6 +108,7 @@ public: static constexpr int tile_size_transform_bi = 4; static constexpr int tile_size_compute_yi = 8; static constexpr int team_size_compute_fused_deidrj = sizeof(real_type) == 4 ? 4 : 2; +#endif // Custom MDRangePolicy, Rank3, to reduce verbosity of kernel launches // This hides the Kokkos::IndexType and Kokkos::Rank<3...> From af231d544702e6edcf1503cde09245530ebbde8a Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 17 Jan 2022 08:54:17 -0700 Subject: [PATCH 17/75] Fix comments in Makefiles --- src/MAKE/MACHINES/Makefile.crusher_kokkos | 2 +- src/MAKE/MACHINES/Makefile.spock_kokkos | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MAKE/MACHINES/Makefile.crusher_kokkos b/src/MAKE/MACHINES/Makefile.crusher_kokkos index 167341ec9c..065c20519f 100644 --- a/src/MAKE/MACHINES/Makefile.crusher_kokkos +++ b/src/MAKE/MACHINES/Makefile.crusher_kokkos @@ -1,4 +1,4 @@ -# summit_kokkos = KOKKOS/CUDA, V100 GPU and Power9, IBM Spectrum MPI, nvcc compiler with gcc +# crusher_kokkos = KOKKOS/HIP, AMD MI250X GPU and AMD EPYC 7A53 “Optimized 3rd Gen EPYC” CPU, Cray MPICH, hipcc compiler SHELL = /bin/sh diff --git a/src/MAKE/MACHINES/Makefile.spock_kokkos b/src/MAKE/MACHINES/Makefile.spock_kokkos index 44fd213962..f1302bf2db 100644 --- a/src/MAKE/MACHINES/Makefile.spock_kokkos +++ b/src/MAKE/MACHINES/Makefile.spock_kokkos @@ -1,4 +1,4 @@ -# summit_kokkos = KOKKOS/CUDA, V100 GPU and Power9, IBM Spectrum MPI, nvcc compiler with gcc +# spock_kokkos = KOKKOS/HIP, AMD MI100 GPU and AMD EPYC 7662 “Rome” CPU, Cray MPICH, hipcc compiler SHELL = /bin/sh From 241a44f1afbed0ccdfb09ede418ea9f0443622ca Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 04:56:47 -0500 Subject: [PATCH 18/75] another workaround for rerun --- src/dump.h | 1 + src/output.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dump.h b/src/dump.h index 35da154d7c..482e87a207 100644 --- a/src/dump.h +++ b/src/dump.h @@ -19,6 +19,7 @@ namespace LAMMPS_NS { class Dump : protected Pointers { + friend class Output; public: char *id; // user-defined name of Dump char *style; // style of Dump diff --git a/src/output.cpp b/src/output.cpp index 449077d1fb..1376d365d8 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -632,7 +632,7 @@ int Output::check_time_dumps(bigint ntimestep) { next_dump_any = MAXBIGINT; for (int idump = 0; idump < ndump; idump++) - if ((last_dump[idump] >= 0) && !update->whichflag) + if ((last_dump[idump] >= 0) && !update->whichflag && !dump[idump]->multifile) error->all(FLERR, "Cannot reset timestep with active dump - must undump first"); if (restart_flag_single) { From f0f4a8e6dc20c1c4afbfa3b758acb50a970fbd2d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 06:01:32 -0500 Subject: [PATCH 19/75] plug some memory leaks in MSM kspace style(s) --- src/KSPACE/msm.cpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 2259d2f829..a980be227d 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -626,15 +626,19 @@ void MSM::allocate() gcall->setup(ngcall_buf1,ngcall_buf2); npergrid = 1; + memory->destroy(gcall_buf1); + memory->destroy(gcall_buf2); memory->create(gcall_buf1,npergrid*ngcall_buf1,"msm:gcall_buf1"); memory->create(gcall_buf2,npergrid*ngcall_buf2,"msm:gcall_buf2"); // allocate memory for each grid level for (int n=0; ndestroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); memory->create3d_offset(qgrid[n],nzlo_out[n],nzhi_out[n], nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:qgrid"); + memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n], nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid"); @@ -654,6 +658,8 @@ void MSM::allocate() gc[n]->setup(ngc_buf1[n],ngc_buf2[n]); npergrid = 1; + memory->destroy(gc_buf1[n]); + memory->destroy(gc_buf2[n]); memory->create(gc_buf1[n],npergrid*ngc_buf1[n],"msm:gc_buf1"); memory->create(gc_buf2[n],npergrid*ngc_buf2[n],"msm:gc_buf2"); } else { @@ -835,11 +841,15 @@ void MSM::deallocate_levels() { if (world_levels) { for (int n=0; n < levels; ++n) { - if (qgrid[n]) + if (qgrid[n]) { memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); + qgrid[n] = nullptr; + } - if (egrid[n]) + if (egrid[n]) { memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); + egrid[n] = nullptr; + } if (gc) { if (gc[n]) { From a02e11040d55e0b99683efd990f374bb446d869d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 06:34:57 -0500 Subject: [PATCH 20/75] whitespace --- src/KSPACE/msm.cpp | 78 +++++++++++++++++++++++----------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index a980be227d..4415175bf0 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -116,7 +116,7 @@ void MSM::settings(int narg, char **arg) MSM::~MSM() { - delete [] factors; + delete[] factors; deallocate(); if (peratom_allocate_flag) deallocate_peratom(); deallocate_levels(); @@ -867,57 +867,57 @@ void MSM::deallocate_levels() } } - delete [] ngrid; + delete[] ngrid; ngrid = nullptr; memory->destroy(procneigh_levels); - delete [] world_levels; - delete [] active_flag; + delete[] world_levels; + delete[] active_flag; - delete [] gc; - delete [] gc_buf1; - delete [] gc_buf2; - delete [] ngc_buf1; - delete [] ngc_buf2; + delete[] gc; + delete[] gc_buf1; + delete[] gc_buf2; + delete[] ngc_buf1; + delete[] ngc_buf2; - delete [] alpha; - delete [] betax; - delete [] betay; - delete [] betaz; + delete[] alpha; + delete[] betax; + delete[] betay; + delete[] betaz; - delete [] nx_msm; - delete [] ny_msm; - delete [] nz_msm; + delete[] nx_msm; + delete[] ny_msm; + delete[] nz_msm; - delete [] nxlo_in; - delete [] nylo_in; - delete [] nzlo_in; + delete[] nxlo_in; + delete[] nylo_in; + delete[] nzlo_in; - delete [] nxhi_in; - delete [] nyhi_in; - delete [] nzhi_in; + delete[] nxhi_in; + delete[] nyhi_in; + delete[] nzhi_in; - delete [] nxlo_out; - delete [] nylo_out; - delete [] nzlo_out; + delete[] nxlo_out; + delete[] nylo_out; + delete[] nzlo_out; - delete [] nxhi_out; - delete [] nyhi_out; - delete [] nzhi_out; + delete[] nxhi_out; + delete[] nyhi_out; + delete[] nzhi_out; - delete [] delxinv; - delete [] delyinv; - delete [] delzinv; + delete[] delxinv; + delete[] delyinv; + delete[] delzinv; - delete [] qgrid; - delete [] egrid; + delete[] qgrid; + delete[] egrid; - delete [] v0grid; - delete [] v1grid; - delete [] v2grid; - delete [] v3grid; - delete [] v4grid; - delete [] v5grid; + delete[] v0grid; + delete[] v1grid; + delete[] v2grid; + delete[] v3grid; + delete[] v4grid; + delete[] v5grid; world_levels = nullptr; active_flag = nullptr; From 752552e0f82b0a979494b9594a5a77e6ae767a7d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 06:39:20 -0500 Subject: [PATCH 21/75] remove duplicate code --- src/KSPACE/msm.cpp | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 4415175bf0..09ea775ef4 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -67,33 +67,7 @@ MSM::MSM(LAMMPS *lmp) factors[0] = 2; MPI_Comm_rank(world,&me); - - phi1d = dphi1d = nullptr; - nmax = 0; - part2grid = nullptr; - - g_direct = nullptr; - g_direct_top = nullptr; - - v0_direct = v1_direct = v2_direct = nullptr; - v3_direct = v4_direct = v5_direct = nullptr; - - v0_direct_top = v1_direct_top = v2_direct_top = nullptr; - v3_direct_top = v4_direct_top = v5_direct_top = nullptr; - - ngrid = nullptr; - - alpha = betax = betay = betaz = nullptr; - nx_msm = ny_msm = nz_msm = nullptr; - nxlo_in = nylo_in = nzlo_in = nullptr; - nxhi_in = nyhi_in = nzhi_in = nullptr; - nxlo_out = nylo_out = nzlo_out = nullptr; - nxhi_out = nyhi_out = nzhi_out = nullptr; - delxinv = delyinv = delzinv = nullptr; - qgrid = nullptr; - egrid = nullptr; - v0grid = v1grid = v2grid = v3grid = v4grid = v5grid = nullptr; peratom_allocate_flag = 0; scalar_pressure_flag = 1; From 2e85233b11936c3b12ff8a8f528967d0a7fb5b5a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 07:20:28 -0500 Subject: [PATCH 22/75] before changing box settings, it must be initialized at least once otherwise change_box would not really be needed, but commands like boundary can (still) be used --- src/change_box.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/change_box.cpp b/src/change_box.cpp index bbac78ab3d..a79b023ac3 100644 --- a/src/change_box.cpp +++ b/src/change_box.cpp @@ -24,6 +24,7 @@ #include "lattice.h" #include "modify.h" #include "output.h" +#include "update.h" #include #include @@ -201,8 +202,11 @@ void ChangeBox::command(int narg, char **arg) scale[0] = domain->lattice->xlattice; scale[1] = domain->lattice->ylattice; scale[2] = domain->lattice->zlattice; - } - else scale[0] = scale[1] = scale[2] = 1.0; + } else scale[0] = scale[1] = scale[2] = 1.0; + + // enforce initialization if there has not one before + + if (!update->first_update) lmp->init(); // perform sequence of operations // first insure atoms are in current box & update box via shrink-wrap @@ -282,9 +286,7 @@ void ChangeBox::command(int narg, char **arg) } else if (ops[m].style == BOUNDARY) { domain->set_boundary(3,&arg[ops[m].boundindex],1); if (domain->dimension == 2 && domain->zperiodic == 0) - error->all(FLERR, - "Cannot change box z boundary to " - "non-periodic for a 2d simulation"); + error->all(FLERR, "Cannot change box z boundary to non-periodic for a 2d simulation"); domain->set_initial_box(); domain->set_global_box(); domain->set_local_box(); From c1f7685c988998a901d1f56bd4237f736d16f518 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 08:30:54 -0500 Subject: [PATCH 23/75] Revert "before changing box settings, it must be initialized at least once" Looking for alternate solution since this change has too many unwanted side effects. This reverts commit 2e85233b11936c3b12ff8a8f528967d0a7fb5b5a. --- src/change_box.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/change_box.cpp b/src/change_box.cpp index a79b023ac3..bbac78ab3d 100644 --- a/src/change_box.cpp +++ b/src/change_box.cpp @@ -24,7 +24,6 @@ #include "lattice.h" #include "modify.h" #include "output.h" -#include "update.h" #include #include @@ -202,11 +201,8 @@ void ChangeBox::command(int narg, char **arg) scale[0] = domain->lattice->xlattice; scale[1] = domain->lattice->ylattice; scale[2] = domain->lattice->zlattice; - } else scale[0] = scale[1] = scale[2] = 1.0; - - // enforce initialization if there has not one before - - if (!update->first_update) lmp->init(); + } + else scale[0] = scale[1] = scale[2] = 1.0; // perform sequence of operations // first insure atoms are in current box & update box via shrink-wrap @@ -286,7 +282,9 @@ void ChangeBox::command(int narg, char **arg) } else if (ops[m].style == BOUNDARY) { domain->set_boundary(3,&arg[ops[m].boundindex],1); if (domain->dimension == 2 && domain->zperiodic == 0) - error->all(FLERR, "Cannot change box z boundary to non-periodic for a 2d simulation"); + error->all(FLERR, + "Cannot change box z boundary to " + "non-periodic for a 2d simulation"); domain->set_initial_box(); domain->set_global_box(); domain->set_local_box(); From 29b5c2659cf8782f79af96f3a684e9cd6cf0bdbe Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 09:04:04 -0500 Subject: [PATCH 24/75] source code formatting --- src/KSPACE/msm.cpp | 6 ++---- src/KSPACE/pppm.cpp | 6 ++---- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 09ea775ef4..9fd581618f 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -1044,8 +1044,7 @@ void MSM::set_grid_global() } if (flag && gridflag && me == 0) - error->warning(FLERR, - "Number of MSM mesh points changed to be a multiple of 2"); + error->warning(FLERR, "Number of MSM mesh points changed to be a multiple of 2"); // adjust Coulombic cutoff to give desired error (if requested) @@ -1071,8 +1070,7 @@ void MSM::set_grid_global() *p_cutoff = cutoff; if (me == 0) - error->warning(FLERR,"Adjusting Coulombic cutoff for " - "MSM, new cutoff = {:.8}", cutoff); + error->warning(FLERR,"Adjusting Coulombic cutoff for MSM, new cutoff = {:.8}", cutoff); } if (triclinic == 0) { diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 06d90e5d0e..d531233a3a 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -197,11 +197,9 @@ void PPPM::init() error->all(FLERR,"Must redefine kspace_style after changing to triclinic box"); if (domain->triclinic && differentiation_flag == 1) - error->all(FLERR,"Cannot (yet) use PPPM with triclinic box " - "and kspace_modify diff ad"); + error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and kspace_modify diff ad"); if (domain->triclinic && slabflag) - error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and " - "slab correction"); + error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and slab correction"); if (domain->dimension == 2) error->all(FLERR,"Cannot use PPPM with 2d simulation"); From 1c8f427e8a4fbcbebb5504aeb1aa095f6337a339 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 09:09:19 -0500 Subject: [PATCH 25/75] detect when MSM::setup() is called before proper initialization and error out --- src/KSPACE/msm.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 9fd581618f..da4a2ad658 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -285,6 +285,11 @@ double MSM::estimate_total_error() void MSM::setup() { + // change_box may trigger MSM::setup() before MSM::init() was called + // error out and request full initialization. + + if (!delxinv) error->all(FLERR, "MSM must be fully initialized for this operation"); + double *prd; double a = cutoff; From 389c35a2d3517573a5b77f96f02dec9ae4aedfb7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 09:09:37 -0500 Subject: [PATCH 26/75] use alternate method to create triclinic box for unit test --- unittest/python/python-commands.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py index bd8512894f..58a0772510 100644 --- a/unittest/python/python-commands.py +++ b/unittest/python/python-commands.py @@ -356,18 +356,16 @@ create_atoms 1 single & def test_extract_box_triclinic(self): self.lmp.command("boundary p p p") - self.lmp.command("region box block 0 2 0 2 0 2") + self.lmp.command("region box prism 0 2 0 2 0 2 0.1 0.2 0.3") self.lmp.command("create_box 1 box") - self.lmp.command("change_box all triclinic") - self.lmp.command("change_box all xy final 0.1 yz final 0.2 xz final 0.3") boxlo, boxhi, xy, yz, xz, periodicity, box_change = self.lmp.extract_box() self.assertEqual(boxlo, [0.0, 0.0, 0.0]) self.assertEqual(boxhi, [2.0, 2.0, 2.0]) self.assertEqual(xy, 0.1) - self.assertEqual(yz, 0.2) - self.assertEqual(xz, 0.3) + self.assertEqual(xz, 0.2) + self.assertEqual(yz, 0.3) self.assertEqual(periodicity, [1, 1, 1]) self.assertEqual(box_change, 0) From 38bcc493b011ca9bc4e1b47efc7fd57870dc73c8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 10:28:40 -0500 Subject: [PATCH 27/75] simplify and modernize CMAP file parser --- src/MOLECULE/fix_cmap.cpp | 148 +++----------------------------------- 1 file changed, 11 insertions(+), 137 deletions(-) diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index a763c5d14c..c6e356e829 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -37,6 +37,7 @@ #include "force.h" #include "math_const.h" #include "memory.h" +#include "potential_file_reader.h" #include "respa.h" #include "tokenizer.h" #include "update.h" @@ -634,150 +635,23 @@ double FixCMAP::compute_scalar() void FixCMAP::read_grid_map(char *cmapfile) { - char linebuf[MAXLINE]; - char *chunk,*line; - int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter; - - FILE *fp = nullptr; if (comm->me == 0) { - fp = utils::open_potential(cmapfile,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open fix cmap file {}: {}", - cmapfile, utils::getsyserror()); + try { + memset(&cmapgrid[0][0][0], 6*CMAPDIM*CMAPDIM, sizeof(double)); + PotentialFileReader reader(lmp, cmapfile, "cmap grid"); - } + // there are six maps in this order. + // alanine, alanine-proline, proline, proline-proline, glycine, glycine-proline. + // read as one big blob of numbers while ignoring comments - for (int ix1 = 0; ix1 < 6; ix1++) - for (int ix2 = 0; ix2 < CMAPDIM; ix2++) - for (int ix3 = 0; ix3 < CMAPDIM; ix3++) - cmapgrid[ix1][ix2][ix3] = 0.0; + reader.next_dvector(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM); - counter = 0; - i1 = i2 = i3 = i4 = i5 = i6 = 0; - j1 = j2 = j3 = j4 = j5 = j6 = 0; - - int done = 0; - - while (!done) { - // only read on rank 0 and broadcast to all other ranks - if (comm->me == 0) - done = (fgets(linebuf,MAXLINE,fp) == nullptr); - - MPI_Bcast(&done,1,MPI_INT,0,world); - if (done) continue; - - MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world); - - // remove leading whitespace - line = linebuf; - while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line; - - // skip if empty line or comment - if (!line || *line =='\n' || *line == '\0' || *line == '#') continue; - - // read in the cmap grid point values - // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors - // will occur if content of the file "cmap.data" is altered - // - // Reading order of the maps: - // 1. Alanine map - // 2. Alanine before proline map - // 3. Proline map - // 4. Two adjacent prolines map - // 5. Glycine map - // 6. Glycine before proline map - - chunk = strtok(line, " \r\n"); - while (chunk != nullptr) { - - // alanine map - - if (counter < CMAPDIM*CMAPDIM) { - cmapgrid[0][i1][j1] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j1++; - if (j1 == CMAPDIM) { - j1 = 0; - i1++; - } - counter++; - } - - // alanine-proline map - - else if (counter >= CMAPDIM*CMAPDIM && - counter < 2*CMAPDIM*CMAPDIM) { - cmapgrid[1][i2][j2]= atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j2++; - if (j2 == CMAPDIM) { - j2 = 0; - i2++; - } - counter++; - } - - // proline map - - else if (counter >= 2*CMAPDIM*CMAPDIM && - counter < 3*CMAPDIM*CMAPDIM) { - cmapgrid[2][i3][j3] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j3++; - if (j3 == CMAPDIM) { - j3 = 0; - i3++; - } - counter++; - } - - // 2 adjacent prolines map - - else if (counter >= 3*CMAPDIM*CMAPDIM && - counter < 4*CMAPDIM*CMAPDIM) { - cmapgrid[3][i4][j4] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j4++; - if (j4 == CMAPDIM) { - j4 = 0; - i4++; - } - counter++; - } - - // glycine map - - else if (counter >= 4*CMAPDIM*CMAPDIM && - counter < 5*CMAPDIM*CMAPDIM) { - cmapgrid[4][i5][j5] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j5++; - if (j5 == CMAPDIM) { - j5 = 0; - i5++; - } - counter++; - } - - // glycine-proline map - - else if (counter >= 5*CMAPDIM*CMAPDIM && - counter < 6*CMAPDIM*CMAPDIM) { - cmapgrid[5][i6][j6] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j6++; - if (j6 == CMAPDIM) { - j6 = 0; - i6++; - } - counter++; - } - - else break; + } catch (std::exception &e) { + error->one(FLERR,"Error reading CMAP potential file: {}", e.what()); } } - if (comm->me == 0) fclose(fp); + MPI_Bcast(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ From be7e6d09399b2d7cbfb13565cb4056e46b71e1fb Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Tue, 18 Jan 2022 12:21:18 -0500 Subject: [PATCH 28/75] Avoid std::string copies --- src/pair_hybrid_scaled.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 5bf593d147..0a4be44be0 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -524,7 +524,7 @@ void PairHybridScaled::write_restart(FILE *fp) int n = scalevars.size(); fwrite(&n, sizeof(int), 1, fp); - for (auto var : scalevars) { + for (auto &var : scalevars) { n = var.size() + 1; fwrite(&n, sizeof(int), 1, fp); fwrite(var.c_str(), sizeof(char), n, fp); From ef13455d6b3ad60cf494b79a588a74a8fb8a4f67 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 16:48:07 -0500 Subject: [PATCH 29/75] add some metadata tags --- examples/PACKAGES/electron_stopping/Si.Si.elstop | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/PACKAGES/electron_stopping/Si.Si.elstop b/examples/PACKAGES/electron_stopping/Si.Si.elstop index 51618149ff..1a6952e571 100644 --- a/examples/PACKAGES/electron_stopping/Si.Si.elstop +++ b/examples/PACKAGES/electron_stopping/Si.Si.elstop @@ -1,8 +1,8 @@ -# Electronic stopping for Si in Si -# Uses metal units +# Electronic stopping for Si in Si UNITS: metal DATE: 2019-02-19 CONTRIBUTOR: Risto Toijala risto.toijala@helsinki.fi # Kinetic energy in eV, stopping power in eV/A # For other atom types, add columns. +# atom type 1 # energy Si in Si 3918.2 6.541 15672.9 13.091 From b890b564caa7571536c750ab16e148018e22fb03 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 16:48:26 -0500 Subject: [PATCH 30/75] modernize electron stopping file reader --- src/EXTRA-FIX/fix_electron_stopping.cpp | 57 +++++++++++-------------- 1 file changed, 26 insertions(+), 31 deletions(-) diff --git a/src/EXTRA-FIX/fix_electron_stopping.cpp b/src/EXTRA-FIX/fix_electron_stopping.cpp index bd54fea97d..49fcdd41bf 100644 --- a/src/EXTRA-FIX/fix_electron_stopping.cpp +++ b/src/EXTRA-FIX/fix_electron_stopping.cpp @@ -29,6 +29,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include "region.h" #include "update.h" @@ -236,49 +237,43 @@ double FixElectronStopping::compute_scalar() return SeLoss_all; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + read electron stopping parameters. only called from MPI rank 0. + format: energy then one column per atom type + read as many lines as available. + energies must be sorted in ascending order. + ---------------------------------------------------------------------- */ void FixElectronStopping::read_table(const char *file) { - char line[MAXLINE]; - - FILE *fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open stopping range table {}: {}", file, utils::getsyserror()); - const int ncol = atom->ntypes + 1; + int nlines = 0; + PotentialFileReader reader(lmp, file, "electron stopping data table"); - int l = 0; - while (true) { - if (fgets(line, MAXLINE, fp) == nullptr) break; // end of file - if (line[0] == '#') continue; // comment + try { + char *line; + double oldvalue = 0.0; - char *pch = strtok(line, " \t\n\r"); - if (pch == nullptr) continue; // blank line + while ((line = reader.next_line())) { + if (nlines >= maxlines) grow_table(); + ValueTokenizer values(line); + elstop_ranges[0][nlines] = values.next_double(); + if (elstop_ranges[0][nlines] <= oldvalue) + throw TokenizerException("energy values must be positive and in ascending order",line); - if (l >= maxlines) grow_table(); + oldvalue = elstop_ranges[0][nlines]; + for (int i = 1; i < ncol; ++i) + elstop_ranges[i][nlines] = values.next_double(); - int i = 0; - for ( ; i < ncol && pch != nullptr; i++) { - elstop_ranges[i][l] = utils::numeric(FLERR, pch,false,lmp); - pch = strtok(nullptr, " \t\n\r"); + ++nlines; } - - if (i != ncol || pch != nullptr) // too short or too long - error->one(FLERR, "fix electron/stopping: Invalid table line"); - - if (l >= 1 && elstop_ranges[0][l] <= elstop_ranges[0][l-1]) - error->one(FLERR, - "fix electron/stopping: Energies must be in ascending order"); - - l++; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing electron stopping data: {}", e.what()); } - table_entries = l; - - if (table_entries == 0) + if (nlines == 0) error->one(FLERR, "Did not find any data in electron/stopping table file"); - fclose(fp); + table_entries = nlines; } /* ---------------------------------------------------------------------- */ From e22f62f6db5ce02e8571150eab62437323e907d4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 19:55:30 -0500 Subject: [PATCH 31/75] add some metadata for fix gle matrix files --- examples/PACKAGES/gle/qt-300k.A | 2 ++ examples/PACKAGES/gle/qt-300k.C | 2 ++ examples/PACKAGES/gle/smart.A | 2 ++ 3 files changed, 6 insertions(+) diff --git a/examples/PACKAGES/gle/qt-300k.A b/examples/PACKAGES/gle/qt-300k.A index 073100ad9d..31ada37b05 100644 --- a/examples/PACKAGES/gle/qt-300k.A +++ b/examples/PACKAGES/gle/qt-300k.A @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 3.904138445158e-4 4.681059722010e-2 3.079778738058e-2 4.428079381336e-2 5.057825203477e-2 2.591193419597e-2 1.513907125942e-2 -4.789343294190e-2 2.037551040972e-2 6.597801861779e-2 -8.528273506602e-3 -2.230839572773e-3 6.593086307870e-3 -6.653653981891e-2 -2.905096373618e-2 -6.597801861779e-2 2.086885297843e-2 1.145471984072e-2 3.111465343867e-2 1.101562087523e-2 -3.264959166486e-2 diff --git a/examples/PACKAGES/gle/qt-300k.C b/examples/PACKAGES/gle/qt-300k.C index 50d68df251..d6239ca8dc 100644 --- a/examples/PACKAGES/gle/qt-300k.C +++ b/examples/PACKAGES/gle/qt-300k.C @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 2.999985914100e+2 5.937593337000e+0 2.144376751500e+2 5.883055908000e+1 -1.119803766000e+2 -6.793381764000e+1 1.379789732400e+1 5.937593337000e+0 3.781851303000e+2 -5.794518522000e+1 -2.178772681500e+2 -1.649310659100e+2 -6.557113452000e+1 3.833830743000e+1 2.144376751500e+2 -5.794518522000e+1 7.325759985000e+2 2.188507713000e+2 -3.704586531000e+2 -3.385193865000e+1 -4.827706758000e+0 diff --git a/examples/PACKAGES/gle/smart.A b/examples/PACKAGES/gle/smart.A index a63bd881c8..8abbdea76a 100644 --- a/examples/PACKAGES/gle/smart.A +++ b/examples/PACKAGES/gle/smart.A @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 4.247358737218e-3 3.231404593065e-2 -9.715629522215e-3 8.199488441225e-3 7.288427565896e-3 -3.634229949055e-3 -3.294395982200e-3 4.887738278128e-2 5.978602802893e-1 -2.079222967202e-1 1.088740438247e-1 4.666954324786e-2 -1.334147134493e-2 -3.914996645811e-2 4.015863091190e-2 2.079222967202e-1 1.645970119444e-1 8.351952991453e-2 5.114740085468e-2 1.217862237677e-2 -4.506711205974e-2 From c4c3b545b2f4c3507dbeb2696c5f8e062a1af037 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 19:55:50 -0500 Subject: [PATCH 32/75] modernize fix gle matrix file reading --- src/EXTRA-FIX/fix_gle.cpp | 104 ++++++++++---------------------------- 1 file changed, 26 insertions(+), 78 deletions(-) diff --git a/src/EXTRA-FIX/fix_gle.cpp b/src/EXTRA-FIX/fix_gle.cpp index de225c6980..d7061e199c 100644 --- a/src/EXTRA-FIX/fix_gle.cpp +++ b/src/EXTRA-FIX/fix_gle.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing authors: Michele Ceriotti (EPFL), Joe Morrone (Stony Brook), - Axel Kohylmeyer (Temple U) + Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ #include "fix_gle.h" @@ -24,6 +24,7 @@ #include "error.h" #include "force.h" #include "memory.h" +#include "potential_file_reader.h" #include "random_mars.h" #include "respa.h" #include "update.h" @@ -221,46 +222,17 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : int seed = utils::inumeric(FLERR,arg[6],false,lmp); // LOADING A matrix - FILE *fgle = nullptr; char *fname = arg[7]; + memset(A, ns1sq, sizeof(double)); if (comm->me == 0) { - fgle = utils::open_potential(fname,lmp,nullptr); - if (fgle == nullptr) - error->one(FLERR,"Cannot open A-matrix file {}: {}",fname, utils::getsyserror()); - utils::logmesg(lmp,"Reading A-matrix from {}\n", fname); - } - - // read each line of the file, skipping blank lines or leading '#' - - char line[MAXLINE],*ptr; - int n,nwords,ndone=0,eof=0; - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fgle); - if (ptr == nullptr) { - eof = 1; - fclose(fgle); - } else n = strlen(line) + 1; + PotentialFileReader reader(lmp,fname,"fix gle A matrix"); + try { + reader.next_dvector(A, ns1sq); + } catch (std::exception &e) { + error->one(FLERR,"Error reading A-matrix: {}", e.what()); } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - - nwords = utils::count_words(line); - if (nwords == 0) continue; - - ptr = strtok(line," \t\n\r\f"); - do { - A[ndone] = atof(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - ndone++; - } while ((ptr != nullptr) && (ndone < ns1sq)); } + MPI_Bcast(A,ns1sq,MPI_DOUBLE,0,world); fnoneq=0; gle_every=1; gle_step=0; for (int iarg=8; iargme == 0) { - fgle = utils::open_potential(fname,lmp,nullptr); - if (fgle == nullptr) - error->one(FLERR,"Cannot open C-matrix file {}: {}",fname, utils::getsyserror()); - utils::logmesg(lmp,"Reading C-matrix from {}\n", fname); - } - - // read each line of the file, skipping blank lines or leading '#' - ndone = eof = 0; - const double cfac = force->boltz / force->mvv2e; - - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fgle); - if (ptr == nullptr) { - eof = 1; - fclose(fgle); - } else n = strlen(line) + 1; + PotentialFileReader reader(lmp,fname,"fix gle C matrix"); + try { + reader.next_dvector(C, ns1sq); + } catch (std::exception &e) { + error->one(FLERR,"Error reading C-matrix: {}", e.what()); } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - - nwords = utils::count_words(line); - if (nwords == 0) continue; - - ptr = strtok(line," \t\n\r\f"); - do { - C[ndone] = cfac*atof(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - ndone++; - } while ((ptr != nullptr) && (ndone < ns1sq)); } + MPI_Bcast(C,ns1sq,MPI_DOUBLE,0,world); + const double cfac = force->boltz / force->mvv2e; + for (int i=0; i < ns1sq; ++i) C[i] *= cfac; } #ifdef GLE_DEBUG @@ -366,12 +314,12 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : FixGLE::~FixGLE() { delete random; - delete [] A; - delete [] C; - delete [] S; - delete [] T; - delete [] TT; - delete [] ST; + delete[] A; + delete[] C; + delete[] S; + delete[] T; + delete[] TT; + delete[] ST; memory->destroy(sqrt_m); memory->destroy(gle_s); From 0595e4d7cd1f676975b96b71f9fa9f36d904cb09 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 00:24:01 -0500 Subject: [PATCH 33/75] refactor parsing of pair style list parameter file --- src/MISC/pair_list.cpp | 258 +++++++++++++++++------------------------ src/MISC/pair_list.h | 20 ++-- 2 files changed, 113 insertions(+), 165 deletions(-) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index d1051ecb11..54afcacaed 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -18,19 +18,28 @@ #include "pair_list.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" -#include "error.h" +#include "text_file_reader.h" +#include "tokenizer.h" +#include +#include +#include +#include using namespace LAMMPS_NS; -static const char * const stylename[] = { - "none", "harmonic", "morse", "lj126", nullptr +enum { NONE = 0, HARM, MORSE, LJ126 }; + +static std::map stylename = { + { "none", NONE }, + { "harmonic", HARM }, + { "morse", MORSE }, + { "lj126", LJ126 } }; // fast power function for integer exponent > 0 @@ -55,7 +64,6 @@ PairList::PairList(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; respa_enable = 0; cut_global = 0.0; - style = nullptr; params = nullptr; check_flag = 1; } @@ -66,7 +74,6 @@ PairList::~PairList() { memory->destroy(setflag); memory->destroy(cutsq); - memory->destroy(style); memory->destroy(params); } @@ -90,7 +97,7 @@ void PairList::compute(int eflag, int vflag) int pc = 0; for (int n=0; n < npairs; ++n) { - const list_parm_t &par = params[n]; + const list_param &par = params[n]; i = atom->map(par.id1); j = atom->map(par.id2); @@ -122,34 +129,34 @@ void PairList::compute(int eflag, int vflag) if (rsq < par.cutsq) { const double r2inv = 1.0/rsq; - if (style[n] == HARM) { + if (par.style == HARM) { const double r = sqrt(rsq); - const double dr = par.parm.harm.r0 - r; - fpair = 2.0*par.parm.harm.k*dr/r; + const double dr = par.param.harm.r0 - r; + fpair = 2.0*par.param.harm.k*dr/r; if (eflag_either) - epair = par.parm.harm.k*dr*dr - par.offset; + epair = par.param.harm.k*dr*dr - par.offset; - } else if (style[n] == MORSE) { + } else if (par.style == MORSE) { const double r = sqrt(rsq); - const double dr = par.parm.morse.r0 - r; - const double dexp = exp(par.parm.morse.alpha * dr); - fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha + const double dr = par.param.morse.r0 - r; + const double dexp = exp(par.param.morse.alpha * dr); + fpair = 2.0*par.param.morse.d0*par.param.morse.alpha * (dexp*dexp - dexp) / r; if (eflag_either) - epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; + epair = par.param.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; - } else if (style[n] == LJ126) { + } else if (par.style == LJ126) { const double r6inv = r2inv*r2inv*r2inv; - const double sig6 = mypow(par.parm.lj126.sigma,6); - fpair = 24.0*par.parm.lj126.epsilon*r6inv + const double sig6 = mypow(par.param.lj126.sigma,6); + fpair = 24.0*par.param.lj126.epsilon*r6inv * (2.0*sig6*sig6*r6inv - sig6) * r2inv; if (eflag_either) - epair = 4.0*par.parm.lj126.epsilon*r6inv + epair = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6) - par.offset; } @@ -210,130 +217,73 @@ void PairList::settings(int narg, char **arg) if (strcmp(arg[2],"check") == 0) check_flag = 1; } - FILE *fp = utils::open_potential(arg[0],lmp,nullptr); - char line[1024]; - if (fp == nullptr) - error->all(FLERR,"Cannot open pair list file"); + std::vector mystyles; + std::vector myparams; - // count lines in file for upper limit of storage needed - int num = 1; - while (fgets(line,1024,fp)) ++num; - rewind(fp); - memory->create(style,num,"pair_list:style"); - memory->create(params,num,"pair_list:params"); - - // now read and parse pair list file for real - npairs = 0; - char *ptr; - tagint id1, id2; - int nharm=0, nmorse=0, nlj126=0; - - while (fgets(line,1024,fp)) { - ptr = strtok(line," \t\n\r\f"); - - // skip empty lines - if (!ptr) continue; - - // skip comment lines starting with # - if (*ptr == '#') continue; - - // get atom ids of pair - id1 = ATOTAGINT(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted pair list file"); - id2 = ATOTAGINT(ptr); - - // get potential type - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted pair list file"); - - style[npairs] = NONE; - list_parm_t &par = params[npairs]; - par.id1 = id1; - par.id2 = id2; - - // harmonic potential - if (strcmp(ptr,stylename[HARM]) == 0) { - style[npairs] = HARM; - - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr == nullptr) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); - par.parm.harm.k = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr == nullptr) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); - par.parm.harm.r0 = utils::numeric(FLERR,ptr,false,lmp); - - ++nharm; - - // morse potential - } else if (strcmp(ptr,stylename[MORSE]) == 0) { - style[npairs] = MORSE; - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.d0 = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.alpha = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.r0 = utils::numeric(FLERR,ptr,false,lmp); - - ++nmorse; - - } else if (strcmp(ptr,stylename[LJ126]) == 0) { - // 12-6 lj potential - style[npairs] = LJ126; - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); - par.parm.lj126.epsilon = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); - par.parm.lj126.sigma = utils::numeric(FLERR,ptr,false,lmp); - - ++nlj126; - - } else { - error->all(FLERR,"Unknown pair list potential style"); - } - - // optional cutoff parameter. if not specified use global value - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr != nullptr) && (*ptr != '#')) { - double cut = utils::numeric(FLERR,ptr,false,lmp); - par.cutsq = cut*cut; - } else { - par.cutsq = cut_global*cut_global; - } - - // count complete entry - ++npairs; - } - fclose(fp); - - // informative output + // read and parse potential file only on MPI rank 0. if (comm->me == 0) { - if (screen) - fprintf(screen,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); - if (logfile) - fprintf(logfile,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); + int nharm, nmorse, nlj126, nskipped; + FILE *fp = utils::open_potential(arg[0],lmp,nullptr); + TextFileReader reader(fp,"pair list coeffs"); + npairs = nharm = nmorse = nlj126 = nskipped = 0; + + try { + char *line; + while ((line = reader.next_line())) { + ValueTokenizer values(line); + list_param oneparam; + oneparam.id1 = values.next_tagint(); + oneparam.id2 = values.next_tagint(); + oneparam.style = stylename[values.next_string()]; + ++npairs; + + switch (oneparam.style) { + + case HARM: + oneparam.param.harm.k = values.next_double(); + oneparam.param.harm.r0 = values.next_double(); + ++nharm; + break; + + case MORSE: + oneparam.param.morse.d0 = values.next_double(); + oneparam.param.morse.alpha = values.next_double(); + oneparam.param.morse.r0 = values.next_double(); + ++nmorse; + break; + + case LJ126: + oneparam.param.lj126.epsilon = values.next_double(); + oneparam.param.lj126.sigma = values.next_double(); + ++nlj126; + break; + + case NONE: // fallthrough + error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}", + utils::trim(line)); + ++nskipped; + break; + } + if (values.has_next()) + oneparam.cutsq = values.next_double(); + else + oneparam.cutsq = cut_global*cut_global; + + myparams.push_back(oneparam); + } + } catch (std::exception &e) { + error->one(FLERR,"Error reading pair list coeffs file: {}", e.what()); + } + utils::logmesg(lmp, "Read {} ({}/{}/{}) interacting pair lines from {}. " + "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, arg[0], nskipped); + + memory->create(params,npairs,"pair_list:params"); + memcpy(params, myparams.data(),npairs*sizeof(list_param)); + fclose(fp); } + MPI_Bcast(&npairs, 1, MPI_INT, 0, world); + if (comm->me != 0) memory->create(params,npairs,"pair_list:params"); + MPI_Bcast(params, npairs*sizeof(list_param), MPI_CHAR, 0, world); } /* ---------------------------------------------------------------------- @@ -374,21 +324,21 @@ void PairList::init_style() if (offset_flag) { for (int n=0; n < npairs; ++n) { - list_parm_t &par = params[n]; + list_param &par = params[n]; - if (style[n] == HARM) { - const double dr = sqrt(par.cutsq) - par.parm.harm.r0; - par.offset = par.parm.harm.k*dr*dr; + if (par.style == HARM) { + const double dr = sqrt(par.cutsq) - par.param.harm.r0; + par.offset = par.param.harm.k*dr*dr; - } else if (style[n] == MORSE) { - const double dr = par.parm.morse.r0 - sqrt(par.cutsq); - const double dexp = exp(par.parm.morse.alpha * dr); - par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp); + } else if (par.style == MORSE) { + const double dr = par.param.morse.r0 - sqrt(par.cutsq); + const double dexp = exp(par.param.morse.alpha * dr); + par.offset = par.param.morse.d0 * (dexp*dexp - 2.0*dexp); - } else if (style[n] == LJ126) { + } else if (par.style == LJ126) { const double r6inv = par.cutsq*par.cutsq*par.cutsq; - const double sig6 = mypow(par.parm.lj126.sigma,6); - par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); + const double sig6 = mypow(par.param.lj126.sigma,6); + par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); } } } @@ -411,7 +361,7 @@ double PairList::init_one(int, int) double PairList::memory_usage() { double bytes = (double)npairs * sizeof(int); - bytes += (double)npairs * sizeof(list_parm_t); + bytes += (double)npairs * sizeof(list_param); const int n = atom->ntypes+1; bytes += (double)n*(n*sizeof(int) + sizeof(int *)); bytes += (double)n*(n*sizeof(double) + sizeof(double *)); diff --git a/src/MISC/pair_list.h b/src/MISC/pair_list.h index 60fd3a7a24..b7161c4f2c 100644 --- a/src/MISC/pair_list.h +++ b/src/MISC/pair_list.h @@ -39,8 +39,6 @@ class PairList : public Pair { protected: void allocate(); - enum { NONE = 0, HARM, MORSE, LJ126 }; - // potential specific parameters struct harm_p { double k, r0; @@ -52,23 +50,23 @@ class PairList : public Pair { double epsilon, sigma; }; - union parm_u { - struct harm_p harm; - struct morse_p morse; - struct lj126_p lj126; + union param_u { + harm_p harm; + morse_p morse; + lj126_p lj126; }; - typedef struct { + struct list_param { + int style; // potential style indicator tagint id1, id2; // global atom ids double cutsq; // cutoff**2 for this pair double offset; // energy offset - union parm_u parm; // parameters for style - } list_parm_t; + union param_u param; // parameters for style + }; protected: double cut_global; // global cutoff distance - int *style; // list of styles for pair interactions - list_parm_t *params; // lisf of pair interaction parameters + list_param *params; // lisf of pair interaction parameters int npairs; // # of atom pairs in global list int check_flag; // 1 if checking for missing pairs }; From f5ad91f9feef5b8179cc5e1dd2e8a12c78579102 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 06:28:53 -0500 Subject: [PATCH 34/75] minor tweaks --- src/MANYBODY/pair_tersoff.cpp | 3 +-- src/MISC/pair_list.cpp | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index f21c3bffca..484adc7436 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -431,8 +431,7 @@ void PairTersoff::read_file(char *file) // transparently convert units for supported conversions int unit_convert = reader.get_unit_convert(); - double conversion_factor = utils::get_conversion_factor(utils::ENERGY, - unit_convert); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert); while ((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index 54afcacaed..5215cdc7e2 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -283,7 +283,7 @@ void PairList::settings(int narg, char **arg) } MPI_Bcast(&npairs, 1, MPI_INT, 0, world); if (comm->me != 0) memory->create(params,npairs,"pair_list:params"); - MPI_Bcast(params, npairs*sizeof(list_param), MPI_CHAR, 0, world); + MPI_Bcast(params, npairs*sizeof(list_param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- From 7676f66674f22296e711a6e358aa93ce424fa09b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 06:29:27 -0500 Subject: [PATCH 35/75] modernize potential file parsing --- src/MANYBODY/pair_edip.cpp | 203 ++++++--------- src/MANYBODY/pair_edip.h | 21 +- src/MANYBODY/pair_edip_multi.cpp | 424 ++++++++++++++----------------- src/MANYBODY/pair_edip_multi.h | 16 +- 4 files changed, 294 insertions(+), 370 deletions(-) diff --git a/src/MANYBODY/pair_edip.cpp b/src/MANYBODY/pair_edip.cpp index af2314b8b9..611b1af455 100644 --- a/src/MANYBODY/pair_edip.cpp +++ b/src/MANYBODY/pair_edip.cpp @@ -31,6 +31,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -368,7 +369,7 @@ void PairEDIP::compute(int eflag, int vflag) directorCos_ik_z = invR_ik * dr_ik[2]; cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + - directorCos_ij_z * directorCos_ik_z; + directorCos_ij_z * directorCos_ik_z; cosTetaDiff = cosTeta + tauFunction; cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; @@ -376,33 +377,33 @@ void PairEDIP::compute(int eflag, int vflag) expMinusQFunctionCosTetaDiffCosTetaDiff = exp(-qFunctionCosTetaDiffCosTetaDiff); potentia3B_factor = lambda * - ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + - eta * qFunctionCosTetaDiffCosTetaDiff); + ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + + eta * qFunctionCosTetaDiffCosTetaDiff); exp3B_ik = preExp3B_ij[neighbor_k]; exp3BDerived_ik = preExp3BDerived_ij[neighbor_k]; forceMod3B_factor1_ij = -exp3BDerived_ij * exp3B_ik * potentia3B_factor; forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * qFunction * cosTetaDiff * - (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); + (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); forceMod3B_factor1_ik = -exp3BDerived_ik * exp3B_ij * potentia3B_factor; forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); forceModCoord += (forceMod3B_factor2 * (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); @@ -765,136 +766,92 @@ double PairEDIP::init_one(int i, int j) void PairEDIP::read_file(char *file) { - int params_per_line = 20; - char **words = new char *[params_per_line + 1]; - memory->sfree(params); params = nullptr; nparams = maxparam = 0; // open file on proc 0 - FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file, lmp, nullptr); - if (fp == nullptr) - error->one(FLERR, "Cannot open EDIP potential file {}: {}", file, utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "edip"); + char *line; - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + while ((line = reader.next_line(NPARAMS_PER_LINE))) { + try { + ValueTokenizer values(line); + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - int n, nwords, ielement, jelement, kelement; - char line[MAXLINE], *ptr; - int eof = 0; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - while (true) { - if (comm->me == 0) { - ptr = fgets(line, MAXLINE, fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else - n = strlen(line) + 1; - } - MPI_Bcast(&eof, 1, MPI_INT, 0, world); - if (eof) break; - MPI_Bcast(&n, 1, MPI_INT, 0, world); - MPI_Bcast(line, n, MPI_CHAR, 0, world); + for (ielement = 0; ielement < nelements; ielement++) + if (iname == elements[ielement]) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (jname == elements[jelement]) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (kname == elements[kelement]) break; + if (kelement == nelements) continue; - // strip comment, skip line if blank + // load up parameter settings and error check their values - if ((ptr = strchr(line, '#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); - // concatenate additional lines until have params_per_line words + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n], MAXLINE - n, fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else - n = strlen(line) + 1; + memset(params + nparams, 0, DELTA*sizeof(Param)); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].A = values.next_double(); + params[nparams].B = values.next_double(); + params[nparams].cutoffA = values.next_double(); + params[nparams].cutoffC = values.next_double(); + params[nparams].alpha = values.next_double(); + params[nparams].beta = values.next_double(); + params[nparams].eta = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lambda = values.next_double(); + params[nparams].mu = values.next_double(); + params[nparams].rho = values.next_double(); + params[nparams].sigma = values.next_double(); + params[nparams].Q0 = values.next_double(); + params[nparams].u1 = values.next_double(); + params[nparams].u2 = values.next_double(); + params[nparams].u3 = values.next_double(); + params[nparams].u4 = values.next_double(); + } catch (std::exception &e) { + error->one(FLERR, "Error reading EDIP potential file: {}", e.what()); } - MPI_Bcast(&eof, 1, MPI_INT, 0, world); - if (eof) break; - MPI_Bcast(&n, 1, MPI_INT, 0, world); - MPI_Bcast(line, n, MPI_CHAR, 0, world); - if ((ptr = strchr(line, '#'))) *ptr = '\0'; - nwords = utils::count_words(line); + + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); + + nparams++; } - - if (nwords != params_per_line) error->all(FLERR, "Incorrect format in EDIP potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line, " \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr, " \t\n\r\f"))) continue; - - // ielement,jelement,kelement = 1st args - // if all 3 args are in element list, then parse this line - // else skip to next entry in file - - for (ielement = 0; ielement < nelements; ielement++) - if (strcmp(words[0], elements[ielement]) == 0) break; - if (ielement == nelements) continue; - for (jelement = 0; jelement < nelements; jelement++) - if (strcmp(words[1], elements[jelement]) == 0) break; - if (jelement == nelements) continue; - for (kelement = 0; kelement < nelements; kelement++) - if (strcmp(words[2], elements[kelement]) == 0) break; - if (kelement == nelements) continue; - - // load up parameter settings and error check their values - - if (nparams == maxparam) { - maxparam += DELTA; - params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params"); - - // make certain all addional allocated storage is initialized - // to avoid false positives when checking with valgrind - - memset(params + nparams, 0, DELTA * sizeof(Param)); - } - - params[nparams].ielement = ielement; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].A = atof(words[3]); - params[nparams].B = atof(words[4]); - params[nparams].cutoffA = atof(words[5]); - params[nparams].cutoffC = atof(words[6]); - params[nparams].alpha = atof(words[7]); - params[nparams].beta = atof(words[8]); - params[nparams].eta = atof(words[9]); - params[nparams].gamm = atof(words[10]); - params[nparams].lambda = atof(words[11]); - params[nparams].mu = atof(words[12]); - params[nparams].rho = atof(words[13]); - params[nparams].sigma = atof(words[14]); - params[nparams].Q0 = atof(words[15]); - params[nparams].u1 = atof(words[16]); - params[nparams].u2 = atof(words[17]); - params[nparams].u3 = atof(words[18]); - params[nparams].u4 = atof(words[19]); - - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || - params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || - params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamm < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || - params[nparams].sigma < 0.0) - error->all(FLERR, "Illegal EDIP parameter"); - - nparams++; } + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete[] words; + if (comm->me != 0) + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -946,7 +903,7 @@ void PairEDIP::setup_params() cutoffC = params[0].cutoffC; sigma = params[0].sigma; lambda = params[0].lambda; - gamm = params[0].gamm; + gamm = params[0].gamma; eta = params[0].eta; Q0 = params[0].Q0; mu = params[0].mu; diff --git a/src/MANYBODY/pair_edip.h b/src/MANYBODY/pair_edip.h index 5812768d55..a0d1e45613 100644 --- a/src/MANYBODY/pair_edip.h +++ b/src/MANYBODY/pair_edip.h @@ -34,14 +34,23 @@ class PairEDIP : public Pair { double init_one(int, int); void init_style(); + static constexpr int NPARAMS_PER_LINE = 20; + protected: struct Param { - double A, B; - double cutoffA, cutoffC, cutsq; - double alpha, beta; - double eta, gamm, lambda, mu, rho, sigma, Q0; - double u1, u2, u3, u4; - int ielement, jelement, kelement; + double A, B; // coefficients for pair interaction I-J + double cutoffA; // cut-off distance for pair interaction I-J + double cutoffC; // lower cut-off distance for calculating Z_I + double alpha; // coefficient for calculating Z_I + double beta; // attractive term for pair I-J + double sigma; // cut-off coefficient for pair I-J + double rho; // pair I-J + double gamma; // coefficient for three-body interaction I-J-K + double eta, lambda; // coefficients for function h(l,Z) + double mu, Q0; // coefficients for function Q(Z) + double u1, u2, u3, u4; // coefficients for function tau(Z) + double cutsq; + int ielement, jelement, kelement, dummy; // dummy added for better alignment }; double *preInvR_ij; diff --git a/src/MANYBODY/pair_edip_multi.cpp b/src/MANYBODY/pair_edip_multi.cpp index 8017fa4f8e..ffbc135a34 100644 --- a/src/MANYBODY/pair_edip_multi.cpp +++ b/src/MANYBODY/pair_edip_multi.cpp @@ -30,6 +30,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -118,8 +119,8 @@ void PairEDIPMulti::compute(int eflag, int vflag) double costheta; double dpairZ,dtripleZ; - // eflag != 0 means compute energy contributions in this step - // vflag != 0 means compute virial contributions in this step + // eflag != 0 means compute energy contributions in this step + // vflag != 0 means compute virial contributions in this step evdwl = 0.0; ev_init(eflag,vflag); @@ -155,39 +156,40 @@ void PairEDIPMulti::compute(int eflag, int vflag) // pre-loop to compute environment coordination f(Z) for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + j = jlist[jj]; + j &= NEIGHMASK; - double delx, dely, delz, r_ij; + double delx, dely, delz, r_ij; - delx = x[j][0] - xtmp; - dely = x[j][1] - ytmp; - delz = x[j][2] - ztmp; - r_ij = delx * delx + dely * dely + delz * delz; + delx = x[j][0] - xtmp; + dely = x[j][1] - ytmp; + delz = x[j][2] - ztmp; + r_ij = delx * delx + dely * dely + delz * delz; - jtype = map[type[j]]; - ijparam = elem3param[itype][jtype][jtype]; - if (r_ij > params[ijparam].cutsq) continue; + jtype = map[type[j]]; + const Param ¶m = params[elem3param[itype][jtype][jtype]]; - r_ij = sqrt(r_ij); + if (r_ij > param.cutsq) continue; - // zeta and its derivative dZ/dr + r_ij = sqrt(r_ij); - if (r_ij < params[ijparam].cutoffC) zeta_i += 1.0; - else { - double f, fdr; - edip_fc(r_ij, ¶ms[ijparam], f, fdr); - zeta_i += f; - dzetair = -fdr / r_ij; + // zeta and its derivative dZ/dr - preForceCoord_counter=numForceCoordPairs*5; - preForceCoord[preForceCoord_counter+0]=dzetair; - preForceCoord[preForceCoord_counter+1]=delx; - preForceCoord[preForceCoord_counter+2]=dely; - preForceCoord[preForceCoord_counter+3]=delz; - preForceCoord[preForceCoord_counter+4]=j; - numForceCoordPairs++; - } + if (r_ij < param.cutoffC) zeta_i += 1.0; + else { + double f, fdr; + edip_fc(r_ij, param, f, fdr); + zeta_i += f; + dzetair = -fdr / r_ij; + + preForceCoord_counter=numForceCoordPairs*5; + preForceCoord[preForceCoord_counter+0]=dzetair; + preForceCoord[preForceCoord_counter+1]=delx; + preForceCoord[preForceCoord_counter+2]=dely; + preForceCoord[preForceCoord_counter+3]=delz; + preForceCoord[preForceCoord_counter+4]=j; + numForceCoordPairs++; + } } // two-body interactions @@ -217,7 +219,7 @@ void PairEDIPMulti::compute(int eflag, int vflag) // already considered in constructing the potential double fdr, fdZ; - edip_pair(r_ij, zeta_i, ¶ms[ijparam], evdwl, fdr, fdZ); + edip_pair(r_ij, zeta_i, params[ijparam], evdwl, fdr, fdZ); fpair = -fdr / r_ij; dpairZ += fdZ; @@ -234,102 +236,102 @@ void PairEDIPMulti::compute(int eflag, int vflag) // three-body Forces for (kk = jj + 1; kk < jnum; kk++) { - double dr_ik[3], r_ik, f_ik[3]; + double dr_ik[3], r_ik, f_ik[3]; - k = jlist[kk]; - k &= NEIGHMASK; - ktype = map[type[k]]; - ikparam = elem3param[itype][ktype][ktype]; - ijkparam = elem3param[itype][jtype][ktype]; + k = jlist[kk]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem3param[itype][ktype][ktype]; + ijkparam = elem3param[itype][jtype][ktype]; - dr_ik[0] = x[k][0] - xtmp; - dr_ik[1] = x[k][1] - ytmp; - dr_ik[2] = x[k][2] - ztmp; - r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + dr_ik[0] = x[k][0] - xtmp; + dr_ik[1] = x[k][1] - ytmp; + dr_ik[2] = x[k][2] - ztmp; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; - if (r_ik > params[ikparam].cutsq) continue; + if (r_ik > params[ikparam].cutsq) continue; - r_ik = sqrt(r_ik); + r_ik = sqrt(r_ik); - costheta=dot3(dr_ij, dr_ik) / r_ij / r_ik; + costheta=dot3(dr_ij, dr_ik) / r_ij / r_ik; - double v1, v2, v3, v4, v5, v6, v7; + double v1, v2, v3, v4, v5, v6, v7; - edip_fcut3(r_ij, ¶ms[ijparam], v1, v2); - edip_fcut3(r_ik, ¶ms[ikparam], v3, v4); - edip_h(costheta, zeta_i, ¶ms[ijkparam], v5, v6, v7); + edip_fcut3(r_ij, params[ijparam], v1, v2); + edip_fcut3(r_ik, params[ikparam], v3, v4); + edip_h(costheta, zeta_i, params[ijkparam], v5, v6, v7); - // potential energy and forces - evdwl = v1 * v3 * v5; - dtripleZ += v1 * v3 * v7; + // potential energy and forces + evdwl = v1 * v3 * v5; + dtripleZ += v1 * v3 * v7; - double dri[3], drj[3], drk[3]; - double dhl, dfr; + double dri[3], drj[3], drk[3]; + double dhl, dfr; - dhl = v1 * v3 * v6; + dhl = v1 * v3 * v6; - costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk); + costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk); - f_ij[0] = -dhl * drj[0]; - f_ij[1] = -dhl * drj[1]; - f_ij[2] = -dhl * drj[2]; - f_ik[0] = -dhl * drk[0]; - f_ik[1] = -dhl * drk[1]; - f_ik[2] = -dhl * drk[2]; + f_ij[0] = -dhl * drj[0]; + f_ij[1] = -dhl * drj[1]; + f_ij[2] = -dhl * drj[2]; + f_ik[0] = -dhl * drk[0]; + f_ik[1] = -dhl * drk[1]; + f_ik[2] = -dhl * drk[2]; - dfr = v2 * v3 * v5; - fpair = -dfr / r_ij; + dfr = v2 * v3 * v5; + fpair = -dfr / r_ij; - f_ij[0] += fpair * dr_ij[0]; - f_ij[1] += fpair * dr_ij[1]; - f_ij[2] += fpair * dr_ij[2]; + f_ij[0] += fpair * dr_ij[0]; + f_ij[1] += fpair * dr_ij[1]; + f_ij[2] += fpair * dr_ij[2]; - dfr = v1 * v4 * v5; - fpair = -dfr / r_ik; + dfr = v1 * v4 * v5; + fpair = -dfr / r_ik; - f_ik[0] += fpair * dr_ik[0]; - f_ik[1] += fpair * dr_ik[1]; - f_ik[2] += fpair * dr_ik[2]; + f_ik[0] += fpair * dr_ik[0]; + f_ik[1] += fpair * dr_ik[1]; + f_ik[2] += fpair * dr_ik[2]; - f[j][0] += f_ij[0]; - f[j][1] += f_ij[1]; - f[j][2] += f_ij[2]; + f[j][0] += f_ij[0]; + f[j][1] += f_ij[1]; + f[j][2] += f_ij[2]; - f[k][0] += f_ik[0]; - f[k][1] += f_ik[1]; - f[k][2] += f_ik[2]; + f[k][0] += f_ik[0]; + f[k][1] += f_ik[1]; + f[k][2] += f_ik[2]; - f[i][0] -= f_ij[0] + f_ik[0]; - f[i][1] -= f_ij[1] + f_ik[1]; - f[i][2] -= f_ij[2] + f_ik[2]; + f[i][0] -= f_ij[0] + f_ik[0]; + f[i][1] -= f_ij[1] + f_ik[1]; + f[i][2] -= f_ij[2] + f_ik[2]; - if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); } } // forces due to environment coordination f(Z) for (int idx = 0; idx < numForceCoordPairs; idx++) { - double delx, dely, delz; + double delx, dely, delz; - preForceCoord_counter = idx * 5; - dzetair = preForceCoord[preForceCoord_counter+0]; - delx = preForceCoord[preForceCoord_counter+1]; - dely = preForceCoord[preForceCoord_counter+2]; - delz = preForceCoord[preForceCoord_counter+3]; - j = static_cast (preForceCoord[preForceCoord_counter+4]); + preForceCoord_counter = idx * 5; + dzetair = preForceCoord[preForceCoord_counter+0]; + delx = preForceCoord[preForceCoord_counter+1]; + dely = preForceCoord[preForceCoord_counter+2]; + delz = preForceCoord[preForceCoord_counter+3]; + j = static_cast (preForceCoord[preForceCoord_counter+4]); - dzetair *= (dpairZ + dtripleZ); + dzetair *= (dpairZ + dtripleZ); - f[j][0] += dzetair * delx; - f[j][1] += dzetair * dely; - f[j][2] += dzetair * delz; + f[j][0] += dzetair * delx; + f[j][1] += dzetair * dely; + f[j][2] += dzetair * delz; - f[i][0] -= dzetair * delx; - f[i][1] -= dzetair * dely; - f[i][2] -= dzetair * delz; + f[i][0] -= dzetair * delx; + f[i][1] -= dzetair * dely; + f[i][2] -= dzetair * delz; - evdwl = 0.0; - if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz); + evdwl = 0.0; + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz); } } @@ -342,13 +344,13 @@ double sqr(double x) } //pair Vij, partial derivatives dVij(r,Z)/dr and dVij(r,Z)/dZ -void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng, +void PairEDIPMulti::edip_pair(double r, double z, const Param ¶m, double &eng, double &fdr, double &fZ) { - double A = param->A; - double B = param->B; - double rho = param->rho; - double beta = param->beta; + double A = param.A; + double B = param.B; + double rho = param.rho; + double beta = param.beta; double v1,v2,v3,v4; v1 = pow(B / r, rho); @@ -361,11 +363,11 @@ void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng, } //function fc(r) in calculating coordination Z and derivative fc'(r) -void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fc(double r, const Param ¶m, double &f, double &fdr) { - double a = param->cutoffA; - double c = param->cutoffC; - double alpha = param->alpha; + double a = param.cutoffA; + double c = param.cutoffC; + double alpha = param.alpha; double x; double v1, v2; @@ -392,10 +394,10 @@ void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr) } //cut-off function for Vij and its derivative fcut2'(r) -void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fcut2(double r, const Param ¶m, double &f, double &fdr) { - double sigma = param->sigma; - double a = param->cutoffA; + double sigma = param.sigma; + double a = param.cutoffA; double v1; if (r > a - 1E-6) @@ -411,12 +413,12 @@ void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr) } //function tau(Z) and its derivative tau'(Z) -void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ) +void PairEDIPMulti::edip_tau(double z, const Param ¶m, double &f, double &fdZ) { - double u1 = param->u1; - double u2 = param->u2; - double u3 = param->u3; - double u4 = param->u4; + double u1 = param.u1; + double u2 = param.u2; + double u3 = param.u3; + double u4 = param.u4; double v1, v2; v1 = exp(-u4 * z); @@ -427,13 +429,13 @@ void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ) } //function h(l,Z) and its partial derivatives dh(l,Z)/dl and dh(l,Z)/dZ -void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f, +void PairEDIPMulti::edip_h(double l, double z, const Param ¶m, double &f, double &fdl, double &fdZ) { - double lambda = param->lambda; - double eta = param->eta; - double Q0 = param->Q0; - double mu = param->mu; + double lambda = param.lambda; + double eta = param.eta; + double Q0 = param.Q0; + double mu = param.mu; double Q, QdZ, Tau, TaudZ; double u2, du2l, du2Z; double v1, v2, v3; @@ -464,10 +466,10 @@ void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f, } //cut-off function for Vijk and its derivative fcut3'(r) -void PairEDIPMulti::edip_fcut3(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fcut3(double r, const Param ¶m, double &f, double &fdr) { - double gamma = param->gamma; - double a = param->cutoffA; + double gamma = param.gamma; + double a = param.cutoffA; double v1; if (r > a - 1E-6) @@ -578,137 +580,92 @@ double PairEDIPMulti::init_one(int i, int j) void PairEDIPMulti::read_file(char *file) { - int params_per_line = 20; - char **words = new char*[params_per_line+1]; - memory->sfree(params); params = nullptr; nparams = maxparam = 0; // open file on proc 0 - FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open EDIP potential file {}: {}",file,utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "edip"); + char *line; - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + while ((line = reader.next_line(NPARAMS_PER_LINE))) { + try { + ValueTokenizer values(line); + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); + for (ielement = 0; ielement < nelements; ielement++) + if (iname == elements[ielement]) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (jname == elements[jelement]) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (kname == elements[kelement]) break; + if (kelement == nelements) continue; - // strip comment, skip line if blank + // load up parameter settings and error check their values - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); - // concatenate additional lines until have params_per_line words + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n],MAXLINE-n,fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; + memset(params + nparams, 0, DELTA*sizeof(Param)); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].A = values.next_double(); + params[nparams].B = values.next_double(); + params[nparams].cutoffA = values.next_double(); + params[nparams].cutoffC = values.next_double(); + params[nparams].alpha = values.next_double(); + params[nparams].beta = values.next_double(); + params[nparams].eta = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lambda = values.next_double(); + params[nparams].mu = values.next_double(); + params[nparams].rho = values.next_double(); + params[nparams].sigma = values.next_double(); + params[nparams].Q0 = values.next_double(); + params[nparams].u1 = values.next_double(); + params[nparams].u2 = values.next_double(); + params[nparams].u3 = values.next_double(); + params[nparams].u4 = values.next_double(); + } catch (std::exception &e) { + error->one(FLERR, "Error reading EDIP potential file: {}", e.what()); } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); + + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); + + nparams++; } - - if (nwords != params_per_line) - error->all(FLERR,"Incorrect format in EDIP potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line," \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue; - - // ielement,jelement,kelement = 1st args - // if all 3 args are in element list, then parse this line - // else skip to next entry in file - - for (ielement = 0; ielement < nelements; ielement++) - if (strcmp(words[0],elements[ielement]) == 0) break; - if (ielement == nelements) continue; - for (jelement = 0; jelement < nelements; jelement++) - if (strcmp(words[1],elements[jelement]) == 0) break; - if (jelement == nelements) continue; - for (kelement = 0; kelement < nelements; kelement++) - if (strcmp(words[2],elements[kelement]) == 0) break; - if (kelement == nelements) continue; - - // load up parameter settings and error check their values - - if (nparams == maxparam) { - maxparam += DELTA; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); - - // make certain all addional allocated storage is initialized - // to avoid false positives when checking with valgrind - - memset(params + nparams, 0, DELTA*sizeof(Param)); - } - - params[nparams].ielement = ielement; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].A = atof(words[3]); - params[nparams].B = atof(words[4]); - params[nparams].cutoffA = atof(words[5]); - params[nparams].cutoffC = atof(words[6]); - params[nparams].alpha = atof(words[7]); - params[nparams].beta = atof(words[8]); - params[nparams].eta = atof(words[9]); - params[nparams].gamma = atof(words[10]); - params[nparams].lambda = atof(words[11]); - params[nparams].mu = atof(words[12]); - params[nparams].rho = atof(words[13]); - params[nparams].sigma = atof(words[14]); - params[nparams].Q0 = atof(words[15]); - params[nparams].u1 = atof(words[16]); - params[nparams].u2 = atof(words[17]); - params[nparams].u3 = atof(words[18]); - params[nparams].u4 = atof(words[19]); - - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || - params[nparams].cutoffA < 0.0 || params[nparams].cutoffC < 0.0 || - params[nparams].alpha < 0.0 || params[nparams].beta < 0.0 || - params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || - params[nparams].rho < 0.0 || params[nparams].sigma < 0.0) - error->all(FLERR,"Illegal EDIP parameter"); - - nparams++; } + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete [] words; + if (comm->me != 0) + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -753,5 +710,4 @@ void PairEDIPMulti::setup() rtmp = sqrt(params[m].cutsq); if (rtmp > cutmax) cutmax = rtmp; } - } diff --git a/src/MANYBODY/pair_edip_multi.h b/src/MANYBODY/pair_edip_multi.h index 3ee7347a56..1070956bc2 100644 --- a/src/MANYBODY/pair_edip_multi.h +++ b/src/MANYBODY/pair_edip_multi.h @@ -34,6 +34,8 @@ class PairEDIPMulti : public Pair { double init_one(int, int); void init_style(); + static constexpr int NPARAMS_PER_LINE = 20; + protected: struct Param { double A, B; // coefficients for pair interaction I-J @@ -48,7 +50,7 @@ class PairEDIPMulti : public Pair { double mu, Q0; // coefficients for function Q(Z) double u1, u2, u3, u4; // coefficients for function tau(Z) double cutsq; - int ielement, jelement, kelement; + int ielement, jelement, kelement, dummy; // dummy added for better alignment }; double *preForceCoord; @@ -63,12 +65,12 @@ class PairEDIPMulti : public Pair { void read_file(char *); void setup(); - void edip_pair(double, double, Param *, double &, double &, double &); - void edip_fc(double, Param *, double &, double &); - void edip_fcut2(double, Param *, double &, double &); - void edip_tau(double, Param *, double &, double &); - void edip_h(double, double, Param *, double &, double &, double &); - void edip_fcut3(double, Param *, double &, double &); + void edip_pair(double, double, const Param &, double &, double &, double &); + void edip_fc(double, const Param &, double &, double &); + void edip_fcut2(double, const Param &, double &, double &); + void edip_tau(double, const Param &, double &, double &); + void edip_h(double, double, const Param &, double &, double &, double &); + void edip_fcut3(double, const Param &, double &, double &); }; } // namespace LAMMPS_NS From c85cdb2732312e07ffaf5acbeeb83f1caffd1772 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 10:11:16 -0500 Subject: [PATCH 36/75] always fall back to using the .so extension if available in the LAMMPS module folder --- python/lammps/core.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/python/lammps/core.py b/python/lammps/core.py index fcd5c76ad5..0bc8d0cb3f 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -122,6 +122,9 @@ class lammps(object): for f in os.listdir(winpath)]): lib_ext = ".dll" modpath = winpath + elif any([f.startswith('liblammps') and f.endswith('.so') + for f in os.listdir(modpath)]): + lib_ext = ".so" else: import platform if platform.system() == "Darwin": From 4f0cde40fd04c8f66cc530056d2a53557708e05a Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Wed, 19 Jan 2022 10:16:53 -0700 Subject: [PATCH 37/75] White no longer exists --- src/MAKE/MACHINES/Makefile.white | 128 ------------------------------- 1 file changed, 128 deletions(-) delete mode 100644 src/MAKE/MACHINES/Makefile.white diff --git a/src/MAKE/MACHINES/Makefile.white b/src/MAKE/MACHINES/Makefile.white deleted file mode 100644 index cb101998b3..0000000000 --- a/src/MAKE/MACHINES/Makefile.white +++ /dev/null @@ -1,128 +0,0 @@ -# white = KOKKOS/CUDA package, OpenMPI with nvcc compiler, Pascal GPU, Power8 CPU - -SHELL = /bin/sh - -# --------------------------------------------------------------------- -# compiler/linker settings -# specify flags and libraries needed for your compiler - -KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) -export OMPI_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper -CC = mpicxx -CCFLAGS = -g -O3 -SHFLAGS = -fPIC -DEPFLAGS = -M - -LINK = mpicxx -LINKFLAGS = -g -O3 -LIB = -SIZE = size - -ARCHIVE = ar -ARFLAGS = -rc -SHLIBFLAGS = -shared -KOKKOS_DEVICES = Cuda -KOKKOS_ARCH = Pascal60,Power8 - -# --------------------------------------------------------------------- -# LAMMPS-specific settings, all OPTIONAL -# specify settings for LAMMPS features you will use -# if you change any -D setting, do full re-compile after "make clean" - -# LAMMPS ifdef settings -# see possible settings in Section 3.5 of the manual - -LMP_INC = -DLAMMPS_GZIP - -# MPI library -# see discussion in Section 3.4 of the manual -# MPI wrapper compiler/linker can provide this info -# can point to dummy MPI library in src/STUBS as in Makefile.serial -# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts -# INC = path for mpi.h, MPI compiler settings -# PATH = path for MPI library -# LIB = name of MPI library - -MPI_INC = -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1 -MPI_PATH = -MPI_LIB = - -# FFT library -# see discussion in Section 3.5.2 of manual -# can be left blank to use provided KISS FFT library -# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings -# PATH = path for FFT library -# LIB = name of FFT library - -FFT_INC = -FFT_PATH = -FFT_LIB = - -# JPEG and/or PNG library -# see discussion in Section 3.5.4 of manual -# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC -# INC = path(s) for jpeglib.h and/or png.h -# PATH = path(s) for JPEG library and/or PNG library -# LIB = name(s) of JPEG library and/or PNG library - -JPG_INC = -JPG_PATH = -JPG_LIB = - -# library for loading shared objects (defaults to -ldl, should be empty on Windows) -# uncomment to change the default - -# override DYN_LIB = - -# --------------------------------------------------------------------- -# build rules and dependencies -# do not edit this section - -include Makefile.package.settings -include Makefile.package - -EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC) -EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH) -EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB) $(DYN_LIB) -EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS) -EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS) - -# Path to src files - -vpath %.cpp .. -vpath %.h .. - -# Link target - -$(EXE): main.o $(LMPLIB) $(EXTRA_LINK_DEPENDS) - $(LINK) $(LINKFLAGS) main.o $(EXTRA_PATH) $(LMPLINK) $(EXTRA_LIB) $(LIB) -o $@ - $(SIZE) $@ - -# Library targets - -$(ARLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) - @rm -f ../$(ARLIB) - $(ARCHIVE) $(ARFLAGS) ../$(ARLIB) $(OBJ) - @rm -f $(ARLIB) - @ln -s ../$(ARLIB) $(ARLIB) - -$(SHLIB): $(OBJ) $(EXTRA_LINK_DEPENDS) - $(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o ../$(SHLIB) \ - $(OBJ) $(EXTRA_LIB) $(LIB) - @rm -f $(SHLIB) - @ln -s ../$(SHLIB) $(SHLIB) - -# Compilation rules - -%.o:%.cpp - $(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $< - -# Individual dependencies - -depend : fastdep.exe $(SRC) - @./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1 - -fastdep.exe: ../DEPEND/fastdep.c - gcc -O -o $@ $< - -sinclude .depend From 0a89e5bd24e76a70401909bba1ad88a0f087bb91 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Wed, 19 Jan 2022 10:18:25 -0700 Subject: [PATCH 38/75] Remove non-ASCII chars --- src/MAKE/MACHINES/Makefile.crusher_kokkos | 2 +- src/MAKE/MACHINES/Makefile.spock_kokkos | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MAKE/MACHINES/Makefile.crusher_kokkos b/src/MAKE/MACHINES/Makefile.crusher_kokkos index 065c20519f..7dc1447d4e 100644 --- a/src/MAKE/MACHINES/Makefile.crusher_kokkos +++ b/src/MAKE/MACHINES/Makefile.crusher_kokkos @@ -1,4 +1,4 @@ -# crusher_kokkos = KOKKOS/HIP, AMD MI250X GPU and AMD EPYC 7A53 “Optimized 3rd Gen EPYC” CPU, Cray MPICH, hipcc compiler +# crusher_kokkos = KOKKOS/HIP, AMD MI250X GPU and AMD EPYC 7A53 "Optimized 3rd Gen EPYC" CPU, Cray MPICH, hipcc compiler SHELL = /bin/sh diff --git a/src/MAKE/MACHINES/Makefile.spock_kokkos b/src/MAKE/MACHINES/Makefile.spock_kokkos index f1302bf2db..a85ebb3039 100644 --- a/src/MAKE/MACHINES/Makefile.spock_kokkos +++ b/src/MAKE/MACHINES/Makefile.spock_kokkos @@ -1,4 +1,4 @@ -# spock_kokkos = KOKKOS/HIP, AMD MI100 GPU and AMD EPYC 7662 “Rome” CPU, Cray MPICH, hipcc compiler +# spock_kokkos = KOKKOS/HIP, AMD MI100 GPU and AMD EPYC 7662 "Rome" CPU, Cray MPICH, hipcc compiler SHELL = /bin/sh From 0a8cbcef2bf7f3261698b273e10868beea9c16cc Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Wed, 19 Jan 2022 10:19:50 -0700 Subject: [PATCH 39/75] Add -DNDEBUG to Kokkos Makefiles --- src/MAKE/MACHINES/Makefile.summit_kokkos | 2 +- src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi | 2 +- src/MAKE/OPTIONS/Makefile.kokkos_mpi_only | 2 +- src/MAKE/OPTIONS/Makefile.kokkos_omp | 2 +- src/MAKE/OPTIONS/Makefile.kokkos_phi | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/MAKE/MACHINES/Makefile.summit_kokkos b/src/MAKE/MACHINES/Makefile.summit_kokkos index 87f8c75da2..f22b27cc74 100644 --- a/src/MAKE/MACHINES/Makefile.summit_kokkos +++ b/src/MAKE/MACHINES/Makefile.summit_kokkos @@ -9,7 +9,7 @@ SHELL = /bin/sh KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) CC = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi b/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi index cb3ef0e442..42a8236c7c 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi +++ b/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi @@ -10,7 +10,7 @@ KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd) export MPICH_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper export OMPI_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only b/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only index 6d5e8d779e..0adb53eef0 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only +++ b/src/MAKE/OPTIONS/Makefile.kokkos_mpi_only @@ -7,7 +7,7 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_omp b/src/MAKE/OPTIONS/Makefile.kokkos_omp index e505da8ae2..82144652dd 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_omp +++ b/src/MAKE/OPTIONS/Makefile.kokkos_omp @@ -7,7 +7,7 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_phi b/src/MAKE/OPTIONS/Makefile.kokkos_phi index b825ad691a..9d5691251c 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_phi +++ b/src/MAKE/OPTIONS/Makefile.kokkos_phi @@ -7,7 +7,7 @@ SHELL = /bin/sh # specify flags and libraries needed for your compiler CC = mpicxx -CCFLAGS = -g -O3 +CCFLAGS = -g -O3 -DNDEBUG SHFLAGS = -fPIC DEPFLAGS = -M From 6e8d9cb532b7ac44fbb9cb461dcc884f136daba4 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Wed, 19 Jan 2022 10:27:02 -0700 Subject: [PATCH 40/75] Port bugfix in Kokkos --- lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp index 6964d5b41b..418aa14c68 100644 --- a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp +++ b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp @@ -694,7 +694,7 @@ std::pair CudaInternal::resize_team_scratch_space( int current_team_scratch = 0; int zero = 0; int one = 1; - while (m_team_scratch_pool[current_team_scratch].compare_exchange_weak( + while (!m_team_scratch_pool[current_team_scratch].compare_exchange_weak( zero, one, std::memory_order_release, std::memory_order_relaxed)) { current_team_scratch = (current_team_scratch + 1) % m_n_team_scratch; } From 0b24c5b498a882a3fcc1aec450f1103f73d2693c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 20 Jan 2022 16:04:07 -0500 Subject: [PATCH 41/75] workaround for segfaults due to NULL pointer dereference with hybrid styles --- src/INTEL/npair_intel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/INTEL/npair_intel.cpp b/src/INTEL/npair_intel.cpp index 395e50006c..76a56ef90f 100644 --- a/src/INTEL/npair_intel.cpp +++ b/src/INTEL/npair_intel.cpp @@ -367,7 +367,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, ty[u] = x[j].y; tz[u] = x[j].z; tjtype[u] = x[j].w; - if (THREE) ttag[u] = tag[j]; + if (THREE && ttag) ttag[u] = tag[j]; } if (FULL == 0 && TRI != 1) { @@ -515,7 +515,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, } if (THREE) { - const tagint jtag = ttag[u]; + const tagint jtag = ttag ? ttag[u] : tag[j]; int flist = 0; if (itag > jtag) { if (((itag+jtag) & 1) == 0) flist = 1; From dbe267cb88ace551797d0258d24cdc183d998606 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 20 Jan 2022 16:04:26 -0500 Subject: [PATCH 42/75] join lines, whitespace --- src/INTEL/fix_intel.cpp | 45 +++++++++++++++++------------------------ 1 file changed, 18 insertions(+), 27 deletions(-) diff --git a/src/INTEL/fix_intel.cpp b/src/INTEL/fix_intel.cpp index 519181be52..2fe6b94515 100644 --- a/src/INTEL/fix_intel.cpp +++ b/src/INTEL/fix_intel.cpp @@ -303,8 +303,7 @@ void FixIntel::init() if (force->pair_match("^hybrid", 0) != nullptr) { _pair_hybrid_flag = 1; if (force->newton_pair != 0 && force->pair->no_virial_fdotr_compute) - error->all(FLERR, - "Intel package requires fdotr virial with newton on."); + error->all(FLERR,"Intel package requires fdotr virial with newton on."); } else _pair_hybrid_flag = 0; @@ -325,8 +324,7 @@ void FixIntel::init() _pair_hybrid_zero = 0; if (force->newton_pair == 0) _pair_hybrid_flag = 0; if (nstyles > 1) - error->all(FLERR, - "Currently, cannot offload more than one intel style with hybrid."); + error->all(FLERR,"Currently, cannot offload more than one intel style with hybrid."); } #endif @@ -356,15 +354,12 @@ void FixIntel::init() void FixIntel::setup(int vflag) { if (neighbor->style != Neighbor::BIN) - error->all(FLERR, - "Currently, neighbor style BIN must be used with Intel package."); + error->all(FLERR,"Currently, neighbor style BIN must be used with Intel package."); if (vflag > 3) - error->all(FLERR, - "Cannot currently get per-atom virials with Intel package."); + error->all(FLERR,"Cannot currently get per-atom virials with Intel package."); #ifdef _LMP_INTEL_OFFLOAD if (neighbor->exclude_setting() != 0) - error->all(FLERR, - "Currently, cannot use neigh_modify exclude with Intel package offload."); + error->all(FLERR,"Currently, cannot use neigh_modify exclude with Intel package offload."); post_force(vflag); #endif } @@ -425,14 +420,14 @@ void FixIntel::pair_init_check(const bool cdmessage) if (_offload_balance != 0.0 && comm->me == 0) { #ifndef __INTEL_COMPILER_BUILD_DATE - error->warning(FLERR, "Unknown Intel Compiler Version\n"); + error->warning(FLERR,"Unknown Intel Compiler Version\n"); #else if (__INTEL_COMPILER_BUILD_DATE != 20131008 && __INTEL_COMPILER_BUILD_DATE < 20141023) - error->warning(FLERR, "Unsupported Intel Compiler."); + error->warning(FLERR,"Unsupported Intel Compiler."); #endif #if !defined(__INTEL_COMPILER) - error->warning(FLERR, "Unsupported Intel Compiler."); + error->warning(FLERR,"Unsupported Intel Compiler."); #endif } @@ -505,8 +500,7 @@ void FixIntel::bond_init_check() { if ((_offload_balance != 0.0) && (atom->molecular != Atom::ATOMIC) && (force->newton_pair != force->newton_bond)) - error->all(FLERR, - "INTEL package requires same setting for newton bond and non-bond."); + error->all(FLERR,"INTEL package requires same setting for newton bond and non-bond."); int intel_pair = 0; if (force->pair_match("/intel$", 0) != nullptr) @@ -517,7 +511,7 @@ void FixIntel::bond_init_check() } if (intel_pair == 0) - error->all(FLERR, "Intel styles for bond/angle/dihedral/improper " + error->all(FLERR,"Intel styles for bond/angle/dihedral/improper " "require intel pair style."); } @@ -534,7 +528,7 @@ void FixIntel::kspace_init_check() } if (intel_pair == 0) - error->all(FLERR, "Intel styles for kspace require intel pair style."); + error->all(FLERR,"Intel styles for kspace require intel pair style."); } /* ---------------------------------------------------------------------- */ @@ -551,7 +545,7 @@ void FixIntel::check_neighbor_intel() _offload_noghost = 0; } if (neighbor->requests[i]->skip && _offload_balance != 0.0) - error->all(FLERR, "Cannot yet use hybrid styles with Intel offload."); + error->all(FLERR,"Cannot yet use hybrid styles with Intel offload."); // avoid flagging a neighbor list as both INTEL and OPENMP if (neighbor->requests[i]->intel) @@ -786,8 +780,7 @@ void FixIntel::add_oresults(const ft * _noalias const f_in, if (f_in[1].w == 1) error->all(FLERR,"Bad matrix inversion in mldivide3"); else - error->all(FLERR, - "Sphere particles not yet supported for gayberne/intel"); + error->all(FLERR,"Sphere particles not yet supported for gayberne/intel"); } } @@ -926,7 +919,7 @@ void FixIntel::add_off_results(const ft * _noalias const f_in, int nlocal = atom->nlocal; if (neighbor->ago == 0) { if (_off_overflow_flag[LMP_OVERFLOW]) - error->one(FLERR, "Neighbor list overflow, boost neigh_modify one"); + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); _offload_nlocal = _off_overflow_flag[LMP_LOCAL_MAX] + 1; _offload_min_ghost = _off_overflow_flag[LMP_GHOST_MIN]; _offload_nghost = _off_overflow_flag[LMP_GHOST_MAX] + 1 - @@ -938,7 +931,7 @@ void FixIntel::add_off_results(const ft * _noalias const f_in, if (atom->torque) if (f_in[1].w < 0.0) - error->all(FLERR, "Bad matrix inversion in mldivide3"); + error->all(FLERR,"Bad matrix inversion in mldivide3"); add_results(f_in, ev_global, _off_results_eatom, _off_results_vatom, 1); // Load balance? @@ -1043,8 +1036,7 @@ void FixIntel::output_timing_data() { timers[TIME_OFFLOAD_PAIR]; double tt = MAX(ht,ct); if (timers[TIME_OFFLOAD_LATENCY] / tt > 0.07 && _separate_coi == 0) - error->warning(FLERR, - "Leaving a core free can improve performance for offload"); + error->warning(FLERR,"Leaving a core free can improve performance for offload"); } fprintf(_tscreen, "------------------------------------------------\n"); } @@ -1109,7 +1101,7 @@ void FixIntel::set_offload_affinity() int ppn = get_ppn(node_rank); if (ppn % _ncops != 0) - error->all(FLERR, "MPI tasks per node must be multiple of offload_cards"); + error->all(FLERR,"MPI tasks per node must be multiple of offload_cards"); ppn = ppn / _ncops; _cop = node_rank / ppn; node_rank = node_rank % ppn; @@ -1214,8 +1206,7 @@ int FixIntel::set_host_affinity(const int nomp) int subscription = nthreads * ppn; if (subscription > ncores) { if (rank == 0) - error->warning(FLERR, - "More MPI tasks/OpenMP threads than available cores"); + error->warning(FLERR,"More MPI tasks/OpenMP threads than available cores"); return 0; } if (subscription == ncores) From 49dca516ea74d8d675f5fb534f790e1550485a57 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 20 Jan 2022 16:48:44 -0500 Subject: [PATCH 43/75] simplify API of Force::store_style() --- src/angle_hybrid.cpp | 2 +- src/bond_hybrid.cpp | 2 +- src/dihedral_hybrid.cpp | 2 +- src/force.cpp | 16 ++++++++-------- src/force.h | 2 +- src/improper_hybrid.cpp | 2 +- src/pair_hybrid.cpp | 2 +- src/pair_hybrid_scaled.cpp | 2 +- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/angle_hybrid.cpp b/src/angle_hybrid.cpp index 12bc9aa97a..41f46ae878 100644 --- a/src/angle_hybrid.cpp +++ b/src/angle_hybrid.cpp @@ -240,7 +240,7 @@ void AngleHybrid::settings(int narg, char **arg) error->all(FLERR, "Angle style hybrid cannot have none as an argument"); styles[nstyles] = force->new_angle(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/bond_hybrid.cpp b/src/bond_hybrid.cpp index dad1010a1e..3139d9ecd8 100644 --- a/src/bond_hybrid.cpp +++ b/src/bond_hybrid.cpp @@ -242,7 +242,7 @@ void BondHybrid::settings(int narg, char **arg) if (strncmp(arg[i], "quartic", 7) == 0) has_quartic = m; styles[nstyles] = force->new_bond(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/dihedral_hybrid.cpp b/src/dihedral_hybrid.cpp index 24a18e6783..0f72f999c9 100644 --- a/src/dihedral_hybrid.cpp +++ b/src/dihedral_hybrid.cpp @@ -241,7 +241,7 @@ void DihedralHybrid::settings(int narg, char **arg) error->all(FLERR, "Dihedral style hybrid cannot have none as an argument"); styles[nstyles] = force->new_dihedral(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/force.cpp b/src/force.cpp index 6b7e9033ca..928b95ee78 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -229,7 +229,7 @@ void Force::create_pair(const std::string &style, int trysuffix) int sflag; pair = new_pair(style,trysuffix,sflag); - store_style(pair_style,style,sflag); + pair_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -350,7 +350,7 @@ void Force::create_bond(const std::string &style, int trysuffix) int sflag; bond = new_bond(style,trysuffix,sflag); - store_style(bond_style,style,sflag); + bond_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -427,7 +427,7 @@ void Force::create_angle(const std::string &style, int trysuffix) int sflag; angle = new_angle(style,trysuffix,sflag); - store_style(angle_style,style,sflag); + angle_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -504,7 +504,7 @@ void Force::create_dihedral(const std::string &style, int trysuffix) int sflag; dihedral = new_dihedral(style,trysuffix,sflag); - store_style(dihedral_style,style,sflag); + dihedral_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -581,7 +581,7 @@ void Force::create_improper(const std::string &style, int trysuffix) int sflag; improper = new_improper(style,trysuffix,sflag); - store_style(improper_style,style,sflag); + improper_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -658,7 +658,7 @@ void Force::create_kspace(const std::string &style, int trysuffix) int sflag; kspace = new_kspace(style,trysuffix,sflag); - store_style(kspace_style,style,sflag); + kspace_style = store_style(style,sflag); } /* ---------------------------------------------------------------------- @@ -729,14 +729,14 @@ KSpace *Force::kspace_match(const std::string &word, int exact) if sflag = 1/2/3, append suffix or suffix2 or suffixp to style ------------------------------------------------------------------------- */ -void Force::store_style(char *&str, const std::string &style, int sflag) +char *Force::store_style(const std::string &style, int sflag) { std::string estyle = style; if (sflag == 1) estyle += std::string("/") + lmp->suffix; else if (sflag == 2) estyle += std::string("/") + lmp->suffix2; else if (sflag == 3) estyle += std::string("/") + lmp->suffixp; - str = utils::strdup(estyle); + return utils::strdup(estyle); } /* ---------------------------------------------------------------------- diff --git a/src/force.h b/src/force.h index 09e13c7bd1..cc9ad06282 100644 --- a/src/force.h +++ b/src/force.h @@ -146,7 +146,7 @@ class Force : protected Pointers { KSpace *new_kspace(const std::string &, int, int &); KSpace *kspace_match(const std::string &, int); - void store_style(char *&, const std::string &, int); + char *store_style(const std::string &, int); void set_special(int, char **); double memory_usage(); diff --git a/src/improper_hybrid.cpp b/src/improper_hybrid.cpp index d3a9403a6b..0354f8e92e 100644 --- a/src/improper_hybrid.cpp +++ b/src/improper_hybrid.cpp @@ -241,7 +241,7 @@ void ImproperHybrid::settings(int narg, char **arg) error->all(FLERR, "Improper style hybrid cannot have none as an argument"); styles[nstyles] = force->new_improper(arg[i], 1, dummy); - force->store_style(keywords[nstyles], arg[i], 0); + keywords[nstyles] = force->store_style(arg[i], 0); istyle = i; if (strcmp(arg[i], "table") == 0) i++; diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index e962e02c9e..8c0326d7b8 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -322,7 +322,7 @@ void PairHybrid::settings(int narg, char **arg) error->all(FLERR,"Pair style hybrid cannot have none as an argument"); styles[nstyles] = force->new_pair(arg[iarg],1,dummy); - force->store_style(keywords[nstyles],arg[iarg],0); + keywords[nstyles] = force->store_style(arg[iarg],0); special_lj[nstyles] = special_coul[nstyles] = nullptr; compute_tally[nstyles] = 1; diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 0a4be44be0..68a6199e19 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -328,7 +328,7 @@ void PairHybridScaled::settings(int narg, char **arg) error->all(FLERR, "Pair style hybrid/scaled cannot have none as an argument"); styles[nstyles] = force->new_pair(arg[iarg], 1, dummy); - force->store_style(keywords[nstyles], arg[iarg], 0); + keywords[nstyles] = force->store_style(arg[iarg], 0); special_lj[nstyles] = special_coul[nstyles] = nullptr; compute_tally[nstyles] = 1; From 4104353d7a2971cf9cde79e20341b26a9a1cbcc6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 20 Jan 2022 16:55:59 -0500 Subject: [PATCH 44/75] plug memory leak --- src/read_data.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/read_data.cpp b/src/read_data.cpp index 89b3ecadec..dda418b8a2 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -903,23 +903,28 @@ void ReadData::command(int narg, char **arg) // restore old styles, when reading with nocoeff flag given if (coeffflag == 0) { - if (force->pair) delete force->pair; + delete force->pair; + delete[] force->pair_style; force->pair = saved_pair; force->pair_style = saved_pair_style; - if (force->bond) delete force->bond; + delete force->bond; + delete[] force->bond_style; force->bond = saved_bond; force->bond_style = saved_bond_style; - if (force->angle) delete force->angle; + delete force->angle; + delete[] force->angle_style; force->angle = saved_angle; force->angle_style = saved_angle_style; - if (force->dihedral) delete force->dihedral; + delete force->dihedral; + delete[] force->dihedral_style; force->dihedral = saved_dihedral; force->dihedral_style = saved_dihedral_style; - if (force->improper) delete force->improper; + delete force->improper; + delete[] force->improper_style; force->improper = saved_improper; force->improper_style = saved_improper_style; From e73158bd7ee77025117ab54b0c79f7043f205ab7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 20 Jan 2022 16:56:12 -0500 Subject: [PATCH 45/75] whitespace/formatting --- src/force.cpp | 28 ++++++++++++++-------------- src/pair_hybrid.cpp | 36 ++++++++++++++++++------------------ 2 files changed, 32 insertions(+), 32 deletions(-) diff --git a/src/force.cpp b/src/force.cpp index 928b95ee78..333a6e4715 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -133,14 +133,14 @@ void _noopt Force::create_factories() Force::~Force() { - delete [] pair_style; - delete [] bond_style; - delete [] angle_style; - delete [] dihedral_style; - delete [] improper_style; - delete [] kspace_style; + delete[] pair_style; + delete[] bond_style; + delete[] angle_style; + delete[] dihedral_style; + delete[] improper_style; + delete[] kspace_style; - delete [] pair_restart; + delete[] pair_restart; if (pair) delete pair; if (bond) delete bond; @@ -220,9 +220,9 @@ void Force::setup() void Force::create_pair(const std::string &style, int trysuffix) { - delete [] pair_style; + delete[] pair_style; if (pair) delete pair; - if (pair_restart) delete [] pair_restart; + if (pair_restart) delete[] pair_restart; pair_style = nullptr; pair = nullptr; pair_restart = nullptr; @@ -345,7 +345,7 @@ char *Force::pair_match_ptr(Pair *ptr) void Force::create_bond(const std::string &style, int trysuffix) { - delete [] bond_style; + delete[] bond_style; if (bond) delete bond; int sflag; @@ -422,7 +422,7 @@ Bond *Force::bond_match(const std::string &style) void Force::create_angle(const std::string &style, int trysuffix) { - delete [] angle_style; + delete[] angle_style; if (angle) delete angle; int sflag; @@ -499,7 +499,7 @@ Angle *Force::angle_match(const std::string &style) void Force::create_dihedral(const std::string &style, int trysuffix) { - delete [] dihedral_style; + delete[] dihedral_style; if (dihedral) delete dihedral; int sflag; @@ -576,7 +576,7 @@ Dihedral *Force::dihedral_match(const std::string &style) void Force::create_improper(const std::string &style, int trysuffix) { - delete [] improper_style; + delete[] improper_style; if (improper) delete improper; int sflag; @@ -653,7 +653,7 @@ Improper *Force::improper_match(const std::string &style) void Force::create_kspace(const std::string &style, int trysuffix) { - delete [] kspace_style; + delete[] kspace_style; if (kspace) delete kspace; int sflag; diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 8c0326d7b8..6287fb97db 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -50,20 +50,20 @@ PairHybrid::~PairHybrid() if (nstyles > 0) { for (int m = 0; m < nstyles; m++) { delete styles[m]; - delete [] keywords[m]; - if (special_lj[m]) delete [] special_lj[m]; - if (special_coul[m]) delete [] special_coul[m]; + delete[] keywords[m]; + if (special_lj[m]) delete[] special_lj[m]; + if (special_coul[m]) delete[] special_coul[m]; } } - delete [] styles; - delete [] keywords; - delete [] multiple; + delete[] styles; + delete[] keywords; + delete[] multiple; - delete [] special_lj; - delete [] special_coul; - delete [] compute_tally; + delete[] special_lj; + delete[] special_coul; + delete[] compute_tally; - delete [] svector; + delete[] svector; if (allocated) { memory->destroy(setflag); @@ -187,7 +187,7 @@ void PairHybrid::compute(int eflag, int vflag) } - delete [] saved_special; + delete[] saved_special; if (vflag_fdotr) virial_fdotr_compute(); } @@ -274,9 +274,9 @@ void PairHybrid::settings(int narg, char **arg) if (nstyles > 0) { for (int m = 0; m < nstyles; m++) { delete styles[m]; - delete [] keywords[m]; - if (special_lj[m]) delete [] special_lj[m]; - if (special_coul[m]) delete [] special_coul[m]; + delete[] keywords[m]; + if (special_lj[m]) delete[] special_lj[m]; + if (special_coul[m]) delete[] special_coul[m]; } delete[] styles; delete[] keywords; @@ -297,8 +297,8 @@ void PairHybrid::settings(int narg, char **arg) // allocate list of sub-styles as big as possibly needed if no extra args - styles = new Pair*[narg]; - keywords = new char*[narg]; + styles = new Pair *[narg]; + keywords = new char *[narg]; multiple = new int[narg]; special_lj = new double*[narg]; @@ -454,7 +454,7 @@ void PairHybrid::init_svector() single_extra = MAX(single_extra,styles[m]->single_extra); if (single_extra) { - delete [] svector; + delete[] svector; svector = new double[single_extra]; } } @@ -667,7 +667,7 @@ void PairHybrid::init_style() neighbor->requests[i]->iskip = iskip; neighbor->requests[i]->ijskip = ijskip; } else { - delete [] iskip; + delete[] iskip; memory->destroy(ijskip); } } From 9c8431bf87dda3d7b87f5acbfc99f7dc79890225 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 06:55:06 -0500 Subject: [PATCH 46/75] correct calls to memset() --- src/EXTRA-FIX/fix_gle.cpp | 5 ++--- src/EXTRA-PAIR/pair_e3b.cpp | 3 +-- src/EXTRA-PAIR/pair_e3b.h | 1 - src/MOLECULE/fix_cmap.cpp | 2 +- src/REAXFF/fix_acks2_reaxff.cpp | 2 +- src/REAXFF/fix_qeq_reaxff.cpp | 2 +- 6 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/EXTRA-FIX/fix_gle.cpp b/src/EXTRA-FIX/fix_gle.cpp index d7061e199c..acb7b84fe1 100644 --- a/src/EXTRA-FIX/fix_gle.cpp +++ b/src/EXTRA-FIX/fix_gle.cpp @@ -211,6 +211,8 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : S = new double[ns1sq]; TT = new double[ns1sq]; ST = new double[ns1sq]; + memset(A,0,sizeof(double)*ns1sq); + memset(C,0,sizeof(double)*ns1sq); // start temperature (t ramp) t_start = utils::numeric(FLERR,arg[4],false,lmp); @@ -223,7 +225,6 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : // LOADING A matrix char *fname = arg[7]; - memset(A, ns1sq, sizeof(double)); if (comm->me == 0) { PotentialFileReader reader(lmp,fname,"fix gle A matrix"); try { @@ -256,14 +257,12 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : if (fnoneq == 0) { t_target=t_start; const double kT = t_target * force->boltz / force->mvv2e; - memset(C,0,sizeof(double)*ns1sq); for (int i=0; ime == 0) { PotentialFileReader reader(lmp,fname,"fix gle C matrix"); try { diff --git a/src/EXTRA-PAIR/pair_e3b.cpp b/src/EXTRA-PAIR/pair_e3b.cpp index 7e865ac6f8..3af8e754c0 100644 --- a/src/EXTRA-PAIR/pair_e3b.cpp +++ b/src/EXTRA-PAIR/pair_e3b.cpp @@ -99,7 +99,7 @@ void PairE3B::compute(int eflag, int vflag) ev_init(eflag, vflag); //clear sumExp array - memset(sumExp, 0.0, nbytes); + memset(sumExp, 0, sizeof(double)*maxID); evdwl = 0.0; pvector[0] = pvector[1] = pvector[2] = pvector[3] = 0.0; @@ -364,7 +364,6 @@ void PairE3B::allocateE3B() maxID = find_maxID(); if (!natoms) error->all(FLERR, "No atoms found"); memory->create(sumExp, maxID, "pair:sumExp"); - nbytes = sizeof(double) * maxID; } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-PAIR/pair_e3b.h b/src/EXTRA-PAIR/pair_e3b.h index f6cf916918..1f38c18293 100644 --- a/src/EXTRA-PAIR/pair_e3b.h +++ b/src/EXTRA-PAIR/pair_e3b.h @@ -50,7 +50,6 @@ class PairE3B : public Pair { int **pairO, ***pairH; // pair lists double ***exps, ****del3, ***fpair3, *sumExp; int maxID; //size of global sumExp array - size_t nbytes; //size of sumExp array in bytes int natoms; //to make sure number of atoms is constant virtual void allocate(); diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index c6e356e829..34f54ef26c 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -637,7 +637,7 @@ void FixCMAP::read_grid_map(char *cmapfile) { if (comm->me == 0) { try { - memset(&cmapgrid[0][0][0], 6*CMAPDIM*CMAPDIM, sizeof(double)); + memset(&cmapgrid[0][0][0], 0, 6*CMAPDIM*CMAPDIM*sizeof(double)); PotentialFileReader reader(lmp, cmapfile, "cmap grid"); // there are six maps in this order. diff --git a/src/REAXFF/fix_acks2_reaxff.cpp b/src/REAXFF/fix_acks2_reaxff.cpp index b6789b1b2e..21cce5b55a 100644 --- a/src/REAXFF/fix_acks2_reaxff.cpp +++ b/src/REAXFF/fix_acks2_reaxff.cpp @@ -429,7 +429,7 @@ void FixACKS2ReaxFF::compute_X() double **x = atom->x; int *mask = atom->mask; - memset(X_diag,0.0,atom->nmax*sizeof(double)); + memset(X_diag,0,atom->nmax*sizeof(double)); // fill in the X matrix m_fill = 0; diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index fa8fa79e00..08a4e21fe5 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -1084,7 +1084,7 @@ void FixQEqReaxFF::vector_add(double* dest, double c, double* v, int k) void FixQEqReaxFF::get_chi_field() { - memset(&chi_field[0],0.0,atom->nmax*sizeof(double)); + memset(&chi_field[0],0,atom->nmax*sizeof(double)); if (!efield) return; const double * const *x = (const double * const *)atom->x; From 67dec40afa7aafccfdb629835e34ad93a58b1b36 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 06:55:20 -0500 Subject: [PATCH 47/75] make certain pointers are initialized --- src/EXTRA-PAIR/pair_harmonic_cut.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-PAIR/pair_harmonic_cut.cpp b/src/EXTRA-PAIR/pair_harmonic_cut.cpp index bd0e48efdc..f2e1c422fe 100644 --- a/src/EXTRA-PAIR/pair_harmonic_cut.cpp +++ b/src/EXTRA-PAIR/pair_harmonic_cut.cpp @@ -35,7 +35,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairHarmonicCut::PairHarmonicCut(LAMMPS *lmp) : Pair(lmp) +PairHarmonicCut::PairHarmonicCut(LAMMPS *lmp) : Pair(lmp), k(nullptr), cut(nullptr) { writedata = 1; } From 3a6bd2e698b06efca393294ac2f212527b4acc1a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 10:47:49 -0500 Subject: [PATCH 48/75] update comment about pair_write restrictions --- doc/src/pair_write.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/pair_write.rst b/doc/src/pair_write.rst index 1d88137255..b62383cb7b 100644 --- a/doc/src/pair_write.rst +++ b/doc/src/pair_write.rst @@ -76,8 +76,10 @@ must be set before this command can be invoked. Due to how the pairwise force is computed, an inner value > 0.0 must be specified even if the potential has a finite value at r = 0.0. -For EAM potentials, the pair_write command only tabulates the -pairwise portion of the potential, not the embedding portion. +The *pair_write* command can only be used for pairwise additive +interactions for which a `Pair::single()` function can be and has +been implemented. This excludes for example manybody potentials +or TIP4P coulomb styles. Related commands """""""""""""""" From 3311b1344cee3676468ebe06679c6b95c28639cb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 10:48:35 -0500 Subject: [PATCH 49/75] consolidate pair-wise vs pairwise spelling --- doc/src/Developer_org.rst | 2 +- doc/src/Errors_warnings.rst | 2 +- doc/src/Run_output.rst | 4 ++-- doc/src/Speed_kokkos.rst | 2 +- doc/src/balance.rst | 2 +- doc/src/compute_pressure_cylinder.rst | 2 +- doc/src/package.rst | 2 +- doc/src/pair_charmm.rst | 2 +- doc/src/pair_dpd.rst | 2 +- doc/src/pair_python.rst | 2 +- doc/src/pair_style.rst | 4 ++-- doc/src/pair_sw.rst | 2 +- doc/src/run_style.rst | 2 +- 13 files changed, 15 insertions(+), 15 deletions(-) diff --git a/doc/src/Developer_org.rst b/doc/src/Developer_org.rst index b5acdf5631..ce36d9a590 100644 --- a/doc/src/Developer_org.rst +++ b/doc/src/Developer_org.rst @@ -225,7 +225,7 @@ follows: commands in an input script. - The Force class computes various forces between atoms. The Pair - parent class is for non-bonded or pair-wise forces, which in LAMMPS + parent class is for non-bonded or pairwise forces, which in LAMMPS also includes many-body forces such as the Tersoff 3-body potential if those are computed by walking pairwise neighbor lists. The Bond, Angle, Dihedral, Improper parent classes are styles for bonded diff --git a/doc/src/Errors_warnings.rst b/doc/src/Errors_warnings.rst index d7a97c1507..2d588a1b77 100644 --- a/doc/src/Errors_warnings.rst +++ b/doc/src/Errors_warnings.rst @@ -416,7 +416,7 @@ This will most likely cause errors in kinetic fluctuations. not defined for the specified atom style. *Molecule has bond topology but no special bond settings* - This means the bonded atoms will not be excluded in pair-wise + This means the bonded atoms will not be excluded in pairwise interactions. *Molecule template for create_atoms has multiple molecules* diff --git a/doc/src/Run_output.rst b/doc/src/Run_output.rst index 8adfd4b293..a988be94ad 100644 --- a/doc/src/Run_output.rst +++ b/doc/src/Run_output.rst @@ -106,7 +106,7 @@ individual ranks. Here is an example output for this section: ---------- The third section above lists the number of owned atoms (Nlocal), -ghost atoms (Nghost), and pair-wise neighbors stored per processor. +ghost atoms (Nghost), and pairwise neighbors stored per processor. The max and min values give the spread of these values across processors with a 10-bin histogram showing the distribution. The total number of histogram counts is equal to the number of processors. @@ -114,7 +114,7 @@ number of histogram counts is equal to the number of processors. ---------- The last section gives aggregate statistics (across all processors) -for pair-wise neighbors and special neighbors that LAMMPS keeps track +for pairwise neighbors and special neighbors that LAMMPS keeps track of (see the :doc:`special_bonds ` command). The number of times neighbor lists were rebuilt is tallied, as is the number of potentially *dangerous* rebuilds. If atom movement triggered neighbor diff --git a/doc/src/Speed_kokkos.rst b/doc/src/Speed_kokkos.rst index 14c2ec680e..8b9b2e99af 100644 --- a/doc/src/Speed_kokkos.rst +++ b/doc/src/Speed_kokkos.rst @@ -214,7 +214,7 @@ threads/task as Nt. The product of these two values should be N, i.e. The default for the :doc:`package kokkos ` command when running on KNL is to use "half" neighbor lists and set the Newton flag to "on" for both pairwise and bonded interactions. This will typically - be best for many-body potentials. For simpler pair-wise potentials, it + be best for many-body potentials. For simpler pairwise potentials, it may be faster to use a "full" neighbor list with Newton flag to "off". Use the "-pk kokkos" :doc:`command-line switch ` to change the default :doc:`package kokkos ` options. See its page for diff --git a/doc/src/balance.rst b/doc/src/balance.rst index 5063a502bb..1d24e467d8 100644 --- a/doc/src/balance.rst +++ b/doc/src/balance.rst @@ -383,7 +383,7 @@ multiple groups, its weight is the product of the weight factors. This weight style is useful in combination with pair style :doc:`hybrid `, e.g. when combining a more costly many-body -potential with a fast pair-wise potential. It is also useful when +potential with a fast pairwise potential. It is also useful when using :doc:`run_style respa ` where some portions of the system have many bonded interactions and others none. It assumes that the computational cost for each group remains constant over time. diff --git a/doc/src/compute_pressure_cylinder.rst b/doc/src/compute_pressure_cylinder.rst index a008254540..9913ef159b 100644 --- a/doc/src/compute_pressure_cylinder.rst +++ b/doc/src/compute_pressure_cylinder.rst @@ -61,7 +61,7 @@ Restrictions This compute currently calculates the pressure tensor contributions for pair styles only (i.e. no bond, angle, dihedral, etc. contributions and in the presence of bonded interactions, the result will be incorrect -due to exclusions for special bonds) and requires pair-wise force +due to exclusions for special bonds) and requires pairwise force calculations not available for most many-body pair styles. K-space calculations are also excluded. Note that this pressure compute outputs the configurational terms only; the kinetic contribution is not included diff --git a/doc/src/package.rst b/doc/src/package.rst index 2cf772ea5a..437601dc60 100644 --- a/doc/src/package.rst +++ b/doc/src/package.rst @@ -460,7 +460,7 @@ using *neigh/thread* *on*, a full neighbor list must also be used. Using is turned on by default only when there are 16K atoms or less owned by an MPI rank and when using a full neighbor list. Not all KOKKOS-enabled potentials support this keyword yet, and only thread over atoms. Many -simple pair-wise potentials such as Lennard-Jones do support threading +simple pairwise potentials such as Lennard-Jones do support threading over both atoms and neighbors. The *newton* keyword sets the Newton flags for pairwise and bonded diff --git a/doc/src/pair_charmm.rst b/doc/src/pair_charmm.rst index d45ef58060..e1469ae323 100644 --- a/doc/src/pair_charmm.rst +++ b/doc/src/pair_charmm.rst @@ -119,7 +119,7 @@ name are the older, original LAMMPS implementations. They compute the LJ and Coulombic interactions with an energy switching function (esw, shown in the formula below as S(r)), which ramps the energy smoothly to zero between the inner and outer cutoff. This can cause -irregularities in pair-wise forces (due to the discontinuous second +irregularities in pairwise forces (due to the discontinuous second derivative of energy at the boundaries of the switching region), which in some cases can result in detectable artifacts in an MD simulation. diff --git a/doc/src/pair_dpd.rst b/doc/src/pair_dpd.rst index 8a8b35e50a..8c61c789a6 100644 --- a/doc/src/pair_dpd.rst +++ b/doc/src/pair_dpd.rst @@ -50,7 +50,7 @@ Style *dpd* computes a force field for dissipative particle dynamics Style *dpd/tstat* invokes a DPD thermostat on pairwise interactions, which is equivalent to the non-conservative portion of the DPD force -field. This pair-wise thermostat can be used in conjunction with any +field. This pairwise thermostat can be used in conjunction with any :doc:`pair style `, and in leiu of per-particle thermostats like :doc:`fix langevin ` or ensemble thermostats like Nose Hoover as implemented by :doc:`fix nvt `. To use diff --git a/doc/src/pair_python.rst b/doc/src/pair_python.rst index 35e07dbd11..65e1cd1611 100644 --- a/doc/src/pair_python.rst +++ b/doc/src/pair_python.rst @@ -164,7 +164,7 @@ Following the *LJCutMelt* example, here are the two functions: .. note:: The evaluation of scripted python code will slow down the - computation pair-wise interactions quite significantly. However, this + computation pairwise interactions quite significantly. However, this can be largely worked around through using the python pair style not for the actual simulation, but to generate tabulated potentials on the fly using the :doc:`pair_write ` command. Please see below diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index f657e29aa3..4f1f1e733f 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -154,10 +154,10 @@ accelerated styles exist. * :doc:`coul/wolf/cs ` - Coulomb via Wolf potential with core/shell adjustments * :doc:`dpd ` - dissipative particle dynamics (DPD) * :doc:`dpd/ext ` - generalized force field for DPD -* :doc:`dpd/ext/tstat ` - pair-wise DPD thermostatting with generalized force field +* :doc:`dpd/ext/tstat ` - pairwise DPD thermostatting with generalized force field * :doc:`dpd/fdt ` - DPD for constant temperature and pressure * :doc:`dpd/fdt/energy ` - DPD for constant energy and enthalpy -* :doc:`dpd/tstat ` - pair-wise DPD thermostatting +* :doc:`dpd/tstat ` - pairwise DPD thermostatting * :doc:`dsmc ` - Direct Simulation Monte Carlo (DSMC) * :doc:`e3b ` - Explicit-three body (E3B) water model * :doc:`drip ` - Dihedral-angle-corrected registry-dependent interlayer potential (DRIP) diff --git a/doc/src/pair_sw.rst b/doc/src/pair_sw.rst index d71999b2d4..d87da43b2c 100644 --- a/doc/src/pair_sw.rst +++ b/doc/src/pair_sw.rst @@ -202,7 +202,7 @@ elements are the same. Thus the two-body parameters for Si interacting with C, comes from the SiCC entry. The three-body parameters can in principle be specific to the three elements of the configuration. In the literature, however, the three-body parameters -are usually defined by simple formulas involving two sets of pair-wise +are usually defined by simple formulas involving two sets of pairwise parameters, corresponding to the ij and ik pairs, where i is the center atom. The user must ensure that the correct combining rule is used to calculate the values of the three-body parameters for diff --git a/doc/src/run_style.rst b/doc/src/run_style.rst index b4d1c22113..fd63c82b90 100644 --- a/doc/src/run_style.rst +++ b/doc/src/run_style.rst @@ -89,7 +89,7 @@ in its 3d FFTs. In this scenario, splitting your P total processors into 2 subsets of processors, P1 in the first partition and P2 in the second partition, can enable your simulation to run faster. This is because the long-range forces in PPPM can be calculated at the same -time as pair-wise and bonded forces are being calculated, and the FFTs +time as pairwise and bonded forces are being calculated, and the FFTs can actually speed up when running on fewer processors. To use this style, you must define 2 partitions where P1 is a multiple From 75d20c40ed4e2c03d5f675f9c201f21688532247 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 11:26:26 -0500 Subject: [PATCH 50/75] remove threebody tag caching altogether since it is not reliable --- src/INTEL/npair_intel.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/INTEL/npair_intel.cpp b/src/INTEL/npair_intel.cpp index 76a56ef90f..9578de9965 100644 --- a/src/INTEL/npair_intel.cpp +++ b/src/INTEL/npair_intel.cpp @@ -310,7 +310,6 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, flt_t * _noalias const tz = ncachez + toffs; int * _noalias const tj = ncachej + toffs; int * _noalias const tjtype = ncachejtype + toffs; - tagint * _noalias const ttag = ncachetag + toffs; flt_t * _noalias itx; flt_t * _noalias ity; @@ -367,7 +366,6 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, ty[u] = x[j].y; tz[u] = x[j].z; tjtype[u] = x[j].w; - if (THREE && ttag) ttag[u] = tag[j]; } if (FULL == 0 && TRI != 1) { @@ -515,7 +513,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, } if (THREE) { - const tagint jtag = ttag ? ttag[u] : tag[j]; + const tagint jtag = tag[j]; int flist = 0; if (itag > jtag) { if (((itag+jtag) & 1) == 0) flist = 1; From 57def1a2d5497c67cc2722437b8f7fc411b07191 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 11:26:55 -0500 Subject: [PATCH 51/75] output reformatting and refactoring --- src/INTEL/fix_intel.cpp | 41 +++++++++++++++++++------------------ src/INTEL/npair_intel.cpp | 4 ++-- src/INTEL/pair_sw_intel.cpp | 3 +-- 3 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/INTEL/fix_intel.cpp b/src/INTEL/fix_intel.cpp index 2fe6b94515..97d4f21788 100644 --- a/src/INTEL/fix_intel.cpp +++ b/src/INTEL/fix_intel.cpp @@ -469,27 +469,28 @@ void FixIntel::pair_init_check(const bool cdmessage) #endif if (_print_pkg_info && comm->me == 0) { - if (screen) { - fprintf(screen, - "----------------------------------------------------------\n"); - if (_offload_balance != 0.0) { - fprintf(screen,"Using Intel Coprocessor with %d threads per core, ", - _offload_tpc); - fprintf(screen,"%d threads per task\n",_offload_threads); - } else { - fprintf(screen,"Using Intel Package without Coprocessor.\n"); - } - fprintf(screen,"Precision: %s\n",kmode); - if (cdmessage) { - #ifdef LMP_USE_AVXCD - fprintf(screen,"AVX512 CD Optimizations: Enabled\n"); - #else - fprintf(screen,"AVX512 CD Optimizations: Disabled\n"); - #endif - } - fprintf(screen, - "----------------------------------------------------------\n"); + utils::logmesg(lmp, "----------------------------------------------------------\n"); + if (_offload_balance != 0.0) { + utils::logmesg(lmp,"Using Intel Coprocessor with {} threads per core, " + "{} threads per task\n",_offload_tpc, _offload_threads); + } else { + utils::logmesg(lmp,"Using INTEL Package without Coprocessor.\n"); } + utils::logmesg(lmp,"Compiler: {}\n",platform::compiler_info()); + #ifdef LMP_SIMD_COMPILER + utils::logmesg(lmp,"SIMD compiler directives: Enabled\n"); + #else + utils::logmesg(lmp,"SIMD compiler directives: Disabled\n"); + #endif + utils::logmesg(lmp,"Precision: {}\n",kmode); + if (cdmessage) { + #ifdef LMP_USE_AVXCD + utils::logmesg(lmp,"AVX512 CD Optimizations: Enabled\n"); + #else + utils::logmesg(lmp,"AVX512 CD Optimizations: Disabled\n"); + #endif + } + utils::logmesg(lmp, "----------------------------------------------------------\n"); } _print_pkg_info = 0; } diff --git a/src/INTEL/npair_intel.cpp b/src/INTEL/npair_intel.cpp index 9578de9965..67cd2517e4 100644 --- a/src/INTEL/npair_intel.cpp +++ b/src/INTEL/npair_intel.cpp @@ -31,8 +31,7 @@ using namespace LAMMPS_NS; NPairIntel::NPairIntel(LAMMPS *lmp) : NPair(lmp) { int ifix = modify->find_fix("package_intel"); if (ifix < 0) - error->all(FLERR, - "The 'package intel' command is required for /intel styles"); + error->all(FLERR,"The 'package intel' command is required for /intel styles"); _fix = static_cast(modify->fix[ifix]); #ifdef _LMP_INTEL_OFFLOAD _cop = _fix->coprocessor_number(); @@ -657,6 +656,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, ns += n2 - pack_offset - maxnbors; #ifdef LMP_INTEL_3BODY_FAST + int alln = n; n = lane; for (int u = pack_offset; u < alln; u++) { neighptr[n] = neighptr2[u]; diff --git a/src/INTEL/pair_sw_intel.cpp b/src/INTEL/pair_sw_intel.cpp index fd043ae5a1..8743fbe54c 100644 --- a/src/INTEL/pair_sw_intel.cpp +++ b/src/INTEL/pair_sw_intel.cpp @@ -1115,8 +1115,7 @@ void PairSWIntel::init_style() int ifix = modify->find_fix("package_intel"); if (ifix < 0) - error->all(FLERR, - "The 'package intel' command is required for /intel styles"); + error->all(FLERR,"The 'package intel' command is required for /intel styles"); fix = static_cast(modify->fix[ifix]); fix->pair_init_check(true); From 6506be94092341f2a114b54ac73b971c5266dca3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Jan 2022 15:56:58 -0500 Subject: [PATCH 52/75] update programming style --- src/INTEL/intel_buffers.cpp | 13 ++++--------- src/fix_balance.cpp | 16 +++++++++------- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/src/INTEL/intel_buffers.cpp b/src/INTEL/intel_buffers.cpp index d86570b0d3..3470c18db8 100644 --- a/src/INTEL/intel_buffers.cpp +++ b/src/INTEL/intel_buffers.cpp @@ -296,9 +296,7 @@ void IntelBuffers::free_list_ptrs() /* ---------------------------------------------------------------------- */ template -void IntelBuffers::grow_data3(NeighList *list, - int *&numneighhalf, - int *&cnumneigh) +void IntelBuffers::grow_data3(NeighList *list, int *&numneighhalf, int *&cnumneigh) { const int size = list->get_maxlocal(); int list_num; @@ -321,10 +319,8 @@ void IntelBuffers::grow_data3(NeighList *list, lmp->memory->destroy(_neigh_list_ptrs[list_num].cnumneigh); lmp->memory->destroy(_neigh_list_ptrs[list_num].numneighhalf); } - lmp->memory->create(_neigh_list_ptrs[list_num].cnumneigh, size, - "_cnumneigh"); - lmp->memory->create(_neigh_list_ptrs[list_num].numneighhalf, size, - "_cnumneigh"); + lmp->memory->create(_neigh_list_ptrs[list_num].cnumneigh, size, "_cnumneigh"); + lmp->memory->create(_neigh_list_ptrs[list_num].numneighhalf, size, "_cnumneigh"); _neigh_list_ptrs[list_num].size = size; } numneighhalf = _neigh_list_ptrs[list_num].numneighhalf; @@ -334,8 +330,7 @@ void IntelBuffers::grow_data3(NeighList *list, /* ---------------------------------------------------------------------- */ template -void IntelBuffers::_grow_list_local(NeighList *list, - const int three_body, +void IntelBuffers::_grow_list_local(NeighList *list, const int three_body, const int offload_end) { free_list_local(); diff --git a/src/fix_balance.cpp b/src/fix_balance.cpp index 2c76f13bf0..6b5c4b6ab1 100644 --- a/src/fix_balance.cpp +++ b/src/fix_balance.cpp @@ -13,20 +13,22 @@ ------------------------------------------------------------------------- */ #include "fix_balance.h" -#include -#include "balance.h" -#include "update.h" + #include "atom.h" +#include "balance.h" #include "comm.h" #include "domain.h" -#include "neighbor.h" -#include "irregular.h" +#include "error.h" +#include "fix_store.h" #include "force.h" +#include "irregular.h" #include "kspace.h" #include "modify.h" -#include "fix_store.h" +#include "neighbor.h" #include "rcb.h" -#include "error.h" +#include "update.h" + +#include using namespace LAMMPS_NS; using namespace FixConst; From 900ff39403a8fa0f24f87fe89e5e075d26a82c99 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 16:33:44 -0500 Subject: [PATCH 53/75] make consistent --- src/tabular_function.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tabular_function.cpp b/src/tabular_function.cpp index a3a904d644..07df49729a 100644 --- a/src/tabular_function.cpp +++ b/src/tabular_function.cpp @@ -1,6 +1,6 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/ Sandia National Laboratories + https://www.lammps.org/, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract From c93fba5e2da9d37079969ab928fa424af94734be Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 16:45:38 -0500 Subject: [PATCH 54/75] make certain that offset is always initialized --- src/MISC/pair_list.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index 5215cdc7e2..f5c82ed5e6 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -232,6 +232,7 @@ void PairList::settings(int narg, char **arg) while ((line = reader.next_line())) { ValueTokenizer values(line); list_param oneparam; + oneparam.offset = 0.0; oneparam.id1 = values.next_tagint(); oneparam.id2 = values.next_tagint(); oneparam.style = stylename[values.next_string()]; From d854e58d002993c6dc627c56bb9801d0e227e0ea Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 16:51:07 -0500 Subject: [PATCH 55/75] add sanity check --- src/reader_native.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/reader_native.cpp b/src/reader_native.cpp index ba7a576a50..065adb5d01 100644 --- a/src/reader_native.cpp +++ b/src/reader_native.cpp @@ -130,6 +130,7 @@ void ReaderNative::skip() // read chunk and skip them read_buf(&nchunk, sizeof(int), 1); + if (nchunk < 0) error->one(FLERR,"Dump file is invalid or corrupted"); int n; for (int i = 0; i < nchunk; i++) { From 458ae4e38aa2df2748c9bba16dbe09238ba2ef4d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 17:06:13 -0500 Subject: [PATCH 56/75] initialize pointers, remove unused class member --- src/EXTRA-MOLECULE/bond_fene_nm.cpp | 2 +- src/EXTRA-MOLECULE/bond_fene_nm.h | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.cpp b/src/EXTRA-MOLECULE/bond_fene_nm.cpp index 147a63512a..8abaf97bec 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.cpp +++ b/src/EXTRA-MOLECULE/bond_fene_nm.cpp @@ -29,7 +29,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondFENENM::BondFENENM(LAMMPS *lmp) : BondFENE(lmp) {} +BondFENENM::BondFENENM(LAMMPS *lmp) : BondFENE(lmp), nn(nullptr), mm(nullptr) {} /* ---------------------------------------------------------------------- */ diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.h b/src/EXTRA-MOLECULE/bond_fene_nm.h index f00394c6d8..4dbf64a331 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.h +++ b/src/EXTRA-MOLECULE/bond_fene_nm.h @@ -38,7 +38,6 @@ class BondFENENM : public BondFENE { virtual void *extract(const char *, int &); protected: - double TWO_1_3; double *nn, *mm; virtual void allocate(); From f814dc5ff9ee7f94391751336443c301c0848641 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 17:06:28 -0500 Subject: [PATCH 57/75] add checks after reading data --- src/reader_native.cpp | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/reader_native.cpp b/src/reader_native.cpp index 065adb5d01..32b2279a60 100644 --- a/src/reader_native.cpp +++ b/src/reader_native.cpp @@ -142,8 +142,7 @@ void ReaderNative::skip() read_lines(2); bigint natoms; int rv = sscanf(line,BIGINT_FORMAT,&natoms); - if (rv != 1) - error->one(FLERR,"Dump file is incorrectly formatted"); + if (rv != 1) error->one(FLERR,"Dump file is incorrectly formatted"); read_lines(5); @@ -164,20 +163,17 @@ void ReaderNative::skip_reading_magic_str() if (is_known_magic_str() && revision > 0x0001) { int len; read_buf(&len, sizeof(int), 1); + if (len < 0) error->one(FLERR,"Dump file is invalid or corrupted"); - if (len > 0) { - // has units - skip_buf(sizeof(char)*len); - } + // has units + if (len > 0) skip_buf(sizeof(char)*len); char flag = 0; read_buf(&flag, sizeof(char), 1); - - if (flag) { - skip_buf(sizeof(double)); - } + if (flag) skip_buf(sizeof(double)); read_buf(&len, sizeof(int), 1); + if (len < 0) error->one(FLERR,"Dump file is invalid or corrupted"); skip_buf(sizeof(char)*len); } } From 81e7583a8d00faba5f2833d8e37ad5cd13e893f8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 17:07:12 -0500 Subject: [PATCH 58/75] initialize pointers, use provided constant, reorder include statements --- src/MOLECULE/bond_quartic.cpp | 25 +++++++++++++------------ src/MOLECULE/bond_quartic.h | 1 - src/OPENMP/bond_quartic_omp.cpp | 8 +++++--- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index b3eea79064..30fd460dcc 100644 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -18,24 +18,25 @@ #include "bond_quartic.h" -#include #include "atom.h" -#include "neighbor.h" #include "comm.h" -#include "force.h" -#include "pair.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "neighbor.h" +#include "pair.h" +#include + +#include "math_const.h" using namespace LAMMPS_NS; +using namespace MathConst; /* ---------------------------------------------------------------------- */ -BondQuartic::BondQuartic(LAMMPS *lmp) : Bond(lmp) -{ - TWO_1_3 = pow(2.0,(1.0/3.0)); -} +BondQuartic::BondQuartic(LAMMPS *lmp) : Bond(lmp), k(nullptr), + b1(nullptr), b2(nullptr), rc(nullptr), u0(nullptr) {} /* ---------------------------------------------------------------------- */ @@ -120,7 +121,7 @@ void BondQuartic::compute(int eflag, int vflag) rb = dr - b2[type]; fbond = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb); - if (rsq < TWO_1_3) { + if (rsq < MY_CUBEROOT2) { sr2 = 1.0/rsq; sr6 = sr2*sr2*sr2; fbond += 48.0*sr6*(sr6-0.5)/rsq; @@ -128,7 +129,7 @@ void BondQuartic::compute(int eflag, int vflag) if (eflag) { ebond = k[type]*r2*ra*rb + u0[type]; - if (rsq < TWO_1_3) ebond += 4.0*sr6*(sr6-1.0) + 1.0; + if (rsq < MY_CUBEROOT2) ebond += 4.0*sr6*(sr6-1.0) + 1.0; } // apply force to each of 2 atoms @@ -336,7 +337,7 @@ double BondQuartic::single(int type, double rsq, int i, int j, eng += k[type]*r2*ra*rb + u0[type]; fforce = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb); - if (rsq < TWO_1_3) { + if (rsq < MY_CUBEROOT2) { sr2 = 1.0/rsq; sr6 = sr2*sr2*sr2; eng += 4.0*sr6*(sr6-1.0) + 1.0; diff --git a/src/MOLECULE/bond_quartic.h b/src/MOLECULE/bond_quartic.h index 35ec705849..3973477e68 100644 --- a/src/MOLECULE/bond_quartic.h +++ b/src/MOLECULE/bond_quartic.h @@ -38,7 +38,6 @@ class BondQuartic : public Bond { double single(int, double, int, int, double &); protected: - double TWO_1_3; double *k, *b1, *b2, *rc, *u0; void allocate(); diff --git a/src/OPENMP/bond_quartic_omp.cpp b/src/OPENMP/bond_quartic_omp.cpp index 00e8357644..427293dd3e 100644 --- a/src/OPENMP/bond_quartic_omp.cpp +++ b/src/OPENMP/bond_quartic_omp.cpp @@ -22,13 +22,15 @@ #include "comm.h" #include "force.h" #include "neighbor.h" - #include "pair.h" #include #include "suffix.h" +#include "math_const.h" + using namespace LAMMPS_NS; +using namespace MathConst; /* ---------------------------------------------------------------------- */ @@ -143,7 +145,7 @@ void BondQuarticOMP::eval(int nfrom, int nto, ThrData * const thr) rb = dr - b2[type]; fbond = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb); - if (rsq < TWO_1_3) { + if (rsq < MY_CUBEROOT2) { sr2 = 1.0/rsq; sr6 = sr2*sr2*sr2; fbond += 48.0*sr6*(sr6-0.5)/rsq; @@ -151,7 +153,7 @@ void BondQuarticOMP::eval(int nfrom, int nto, ThrData * const thr) if (EFLAG) { ebond = k[type]*r2*ra*rb + u0[type]; - if (rsq < TWO_1_3) ebond += 4.0*sr6*(sr6-1.0) + 1.0; + if (rsq < MY_CUBEROOT2) ebond += 4.0*sr6*(sr6-1.0) + 1.0; } // apply force to each of 2 atoms From cd165562569770015214004cf49b320c3fe423f7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Jan 2022 17:13:54 -0500 Subject: [PATCH 59/75] add missing break statements. @GenieTim this bug may have tainted your results. you would always get the energy value for any of the dx, dy, dz keywords --- src/compute_pair_local.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index ff9acdc4ef..b98f6eec33 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -277,10 +277,13 @@ int ComputePairLocal::compute_pairs(int flag) break; case DX: ptr[n] = delx*directionCorrection; + break; case DY: ptr[n] = dely*directionCorrection; + break; case DZ: ptr[n] = delz*directionCorrection; + break; case ENG: ptr[n] = eng; break; From 88dadd0fffea798fcd56a9253cb6d020ea41f83c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Jan 2022 21:12:22 -0500 Subject: [PATCH 60/75] consistently refer to INTEL package in upper case --- src/INTEL/fix_intel.cpp | 13 ++++++------- src/INTEL/nbin_intel.cpp | 6 ++---- src/INTEL/verlet_lrt_intel.cpp | 2 +- 3 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/INTEL/fix_intel.cpp b/src/INTEL/fix_intel.cpp index 97d4f21788..6a848a366e 100644 --- a/src/INTEL/fix_intel.cpp +++ b/src/INTEL/fix_intel.cpp @@ -43,7 +43,7 @@ using namespace FixConst; #ifdef __INTEL_OFFLOAD #ifndef _LMP_INTEL_OFFLOAD -#warning "Not building Intel package with Xeon Phi offload support." +#warning "Not building INTEL package with Xeon Phi offload support." #endif #endif @@ -303,7 +303,7 @@ void FixIntel::init() if (force->pair_match("^hybrid", 0) != nullptr) { _pair_hybrid_flag = 1; if (force->newton_pair != 0 && force->pair->no_virial_fdotr_compute) - error->all(FLERR,"Intel package requires fdotr virial with newton on."); + error->all(FLERR,"INTEL package requires fdotr virial with newton on."); } else _pair_hybrid_flag = 0; @@ -354,12 +354,12 @@ void FixIntel::init() void FixIntel::setup(int vflag) { if (neighbor->style != Neighbor::BIN) - error->all(FLERR,"Currently, neighbor style BIN must be used with Intel package."); + error->all(FLERR,"Currently, neighbor style BIN must be used with INTEL package."); if (vflag > 3) - error->all(FLERR,"Cannot currently get per-atom virials with Intel package."); + error->all(FLERR,"Cannot currently get per-atom virials with INTEL package."); #ifdef _LMP_INTEL_OFFLOAD if (neighbor->exclude_setting() != 0) - error->all(FLERR,"Currently, cannot use neigh_modify exclude with Intel package offload."); + error->all(FLERR,"Currently, cannot use neigh_modify exclude with INTEL package offload."); post_force(vflag); #endif } @@ -512,8 +512,7 @@ void FixIntel::bond_init_check() } if (intel_pair == 0) - error->all(FLERR,"Intel styles for bond/angle/dihedral/improper " - "require intel pair style."); + error->all(FLERR,"Intel styles for bond/angle/dihedral/improper require intel pair style."); } /* ---------------------------------------------------------------------- */ diff --git a/src/INTEL/nbin_intel.cpp b/src/INTEL/nbin_intel.cpp index 94f18002a0..29bc7c9ced 100644 --- a/src/INTEL/nbin_intel.cpp +++ b/src/INTEL/nbin_intel.cpp @@ -161,10 +161,8 @@ void NBinIntel::bin_atoms(IntelBuffers * buffers) { const flt_t dx = (INTEL_BIGP - bboxhi[0]); const flt_t dy = (INTEL_BIGP - bboxhi[1]); const flt_t dz = (INTEL_BIGP - bboxhi[2]); - if (dx * dx + dy * dy + dz * dz < - static_cast(neighbor->cutneighmaxsq)) - error->one(FLERR, - "Intel package expects no atoms within cutoff of {1e15,1e15,1e15}."); + if (dx * dx + dy * dy + dz * dz < static_cast(neighbor->cutneighmaxsq)) + error->one(FLERR,"INTEL package expects no atoms within cutoff of (1e15,1e15,1e15)."); } // ---------- Grow and cast/pack buffers ------------- diff --git a/src/INTEL/verlet_lrt_intel.cpp b/src/INTEL/verlet_lrt_intel.cpp index 216ba98302..f867ad73e6 100644 --- a/src/INTEL/verlet_lrt_intel.cpp +++ b/src/INTEL/verlet_lrt_intel.cpp @@ -71,7 +71,7 @@ void VerletLRTIntel::init() #ifndef LMP_INTEL_USELRT error->all(FLERR, - "LRT otion for Intel package disabled at compile time"); + "LRT otion for INTEL package disabled at compile time"); #endif } From f9a2006d732670265cbb0d47631bb3d4291150b3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 24 Jan 2022 08:44:29 -0500 Subject: [PATCH 61/75] reorder includes --- src/ASPHERE/pair_gayberne.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp index 08af6ee3c9..39cb88f81b 100644 --- a/src/ASPHERE/pair_gayberne.cpp +++ b/src/ASPHERE/pair_gayberne.cpp @@ -18,18 +18,18 @@ #include "pair_gayberne.h" -#include -#include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" -#include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "citeme.h" -#include "memory.h" +#include "comm.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include using namespace LAMMPS_NS; From 16810b84eba91c85e23dd9d0f62392d0d7f4ec87 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 14:35:27 -0700 Subject: [PATCH 62/75] Add option to balance dump file output --- src/dump.cpp | 165 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/dump.h | 2 + 2 files changed, 167 insertions(+) diff --git a/src/dump.cpp b/src/dump.cpp index e5652c210c..cdf2f57d9e 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -75,6 +75,7 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp) clearstep = 0; sort_flag = 0; + balance_flag = 0; append_flag = 0; buffer_allow = 0; buffer_flag = 0; @@ -222,6 +223,9 @@ void Dump::init() ids = idsort = nullptr; index = proclist = nullptr; irregular = nullptr; + + if (balance_flag) + error->all(FLERR,"Cannot balance dump output without sorting enabled"); } if (sort_flag) { @@ -417,6 +421,7 @@ void Dump::write() if (sort_flag && sortcol == 0) pack(ids); else pack(nullptr); if (sort_flag) sort(); + if (balance_flag) balance(); // write timestep header // for multiproc, @@ -888,6 +893,161 @@ int Dump::bufcompare_reverse(const int i, const int j, void *ptr) #endif +/* ---------------------------------------------------------------------- + parallel load balance of buf across all procs + must come after sort +------------------------------------------------------------------------- */ + +void Dump::balance() +{ + bigint *proc_offsets,*proc_new_offsets; + memory->create(proc_offsets,nprocs+1,"dump:proc_offsets"); + memory->create(proc_new_offsets,nprocs+1,"dump:proc_new_offsets"); + + // compute atom offset for this proc + + bigint offset; + bigint bnme = nme; + MPI_Scan(&bnme,&offset,1,MPI_LMP_BIGINT,MPI_SUM,world); + + // gather atom offsets for all procs + + MPI_Allgather(&offset,1,MPI_LMP_BIGINT,&proc_offsets[1],1,MPI_LMP_BIGINT,world); + + proc_offsets[0] = 0; + + // how many atoms should I own after balance + + int nme_balance = static_cast(ntotal/nprocs); + + // include remainder atoms on first procs + + int remainder = ntotal % nprocs; + if (me < remainder) nme_balance += 1; + + // compute new atom offset for this proc + + bigint offset_balance; + bigint bnme_balance = nme_balance; + MPI_Scan(&bnme_balance,&offset_balance,1,MPI_LMP_BIGINT,MPI_SUM,world); + + // gather new atom offsets for all procs + + MPI_Allgather(&offset_balance,1,MPI_LMP_BIGINT,&proc_new_offsets[1],1,MPI_LMP_BIGINT,world); + + proc_new_offsets[0] = 0; + + bigint start = proc_new_offsets[me]; + bigint end = proc_new_offsets[me+1]; + + // reset buf size to largest of any post-balance nme values + // this insures proc 0 can receive everyone's info + + int nmax; + MPI_Allreduce(&nme_balance,&nmax,1,MPI_INT,MPI_MAX,world); + + // allocate a second buffer for balanced data + + double* buf_balance; + memory->create(buf_balance,nmax*size_one,"dump:buf_balance"); + + // compute from which procs I am receiving atoms + // post recvs first + + int nswap = 0; + MPI_Request *request = new MPI_Request[nprocs]; + int procstart = start; + int iproc = me; + int iproc_prev; + + for (bigint i = start; i < end; i++) { + + // find which proc this atom belongs to + + while (i < proc_offsets[iproc]) iproc--; + while (i > proc_offsets[iproc+1]-1) iproc++; + + if (i != start && (iproc != iproc_prev || i == end-1)) { + + // finished with proc + + int procrecv = iproc; + if (iproc != iproc_prev) procrecv = iproc_prev; + + int procnrecv = i - procstart + 1; + if (iproc != iproc_prev) procnrecv--; + + // post receive for this proc + + if (iproc_prev != me) + MPI_Irecv(&buf_balance[(procstart - start)*size_one],procnrecv*size_one,MPI_DOUBLE, + procrecv,0,world,&request[nswap++]); + + procstart = i; + } + + iproc_prev = iproc; + } + + // compute which atoms I am sending and to which procs + + procstart = 0; + iproc = me; + for (int i = 0; i < nme; i++) { + + // find which proc this atom should belong to + + while (proc_offsets[me] + i < proc_new_offsets[iproc]) iproc--; + while (proc_offsets[me] + i > proc_new_offsets[iproc+1] - 1) iproc++; + + if (i != 0 && (iproc != iproc_prev || i == nme - 1)) { + + // finished with proc + + int procsend = iproc; + if (iproc != iproc_prev) procsend = iproc_prev; + + int procnsend = i - procstart + 1; + if (iproc != iproc_prev) procnsend--; + + // send for this proc + + if (iproc_prev != me) { + MPI_Send(&buf[procstart*size_one],procnsend*size_one,MPI_DOUBLE,procsend,0,world); + } else { + + // sending to self, copy buffers + + int offset_me = proc_offsets[me] - start; + memcpy(&buf_balance[(offset_me + procstart)*size_one],&buf[procstart*size_one],procnsend*size_one*sizeof(double)); + } + + procstart = i; + } + + iproc_prev = iproc; + } + + // wait for all recvs + + for (int n = 0; n < nswap; n++) + MPI_Wait(&request[n],MPI_STATUS_IGNORE); + + nme = nme_balance; + + // swap buffers + + double *tmp = buf; + buf = buf_balance; + + // cleanup + + memory->destroy(tmp); + memory->destroy(proc_offsets); + memory->destroy(proc_new_offsets); + delete [] request; +} + /* ---------------------------------------------------------------------- process params common to all dumps here if unknown param, call modify_param specific to the dump @@ -1108,6 +1268,11 @@ void Dump::modify_params(int narg, char **arg) } iarg += 2; + } else if (strcmp(arg[iarg],"balance") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); + balance_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + iarg += 2; + } else if (strcmp(arg[iarg],"time") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); time_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); diff --git a/src/dump.h b/src/dump.h index 482e87a207..83edda3810 100644 --- a/src/dump.h +++ b/src/dump.h @@ -67,6 +67,7 @@ class Dump : protected Pointers { int header_flag; // 0 = item, 2 = xyz int flush_flag; // 0 if no flush, 1 if flush every dump int sort_flag; // 1 if sorted output + int balance_flag; // 1 if load balanced output int append_flag; // 1 if open file in append mode, 0 if not int buffer_allow; // 1 if style allows for buffer_flag, 0 if not int buffer_flag; // 1 if buffer output as one big string, 0 if not @@ -161,6 +162,7 @@ class Dump : protected Pointers { static int bufcompare(const int, const int, void *); static int bufcompare_reverse(const int, const int, void *); #endif + void balance(); }; } // namespace LAMMPS_NS From 1d60e4c4639614cfc95e1e291e030d1398fc16da Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 15:16:49 -0700 Subject: [PATCH 63/75] Add docs --- doc/src/dump_modify.rst | 11 ++++++++++- src/dump.cpp | 3 ++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index be75153f6f..52efa47d7c 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -17,13 +17,14 @@ Syntax * one or more keyword/value pairs may be appended * these keywords apply to various dump styles -* keyword = *append* or *at* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* +* keyword = *append* or *at* or *balance* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* .. parsed-literal:: *append* arg = *yes* or *no* *at* arg = N N = index of frame written upon first dump + *balance* arg = *yes* or *no* *buffer* arg = *yes* or *no* *delay* arg = Dstep Dstep = delay output until this timestep @@ -658,6 +659,13 @@ atom ID. A sort value of N or -N means sort the output by the value in the Nth column of per-atom info in either ascending or descending order. +In a parallel run, the per-processor dump file pieces can have +significant imbalance in number of lines of per-atom info. The *balance* +keyword determines whether lines of per-atom output for each processor +snapshot in a parallel run are balanced to be nearly the same. A balance +value of *no* means no balancing will be done, while *yes* means +balancing will be performed. For a serial run, this option is ignored. + The dump *local* style cannot be sorted by atom ID, since there are typically multiple lines of output per atom. Some dump styles, such as *dcd* and *xtc*, require sorting by atom ID to format the output @@ -832,6 +840,7 @@ Default The option defaults are * append = no +* balance = no * buffer = yes for dump styles *atom*, *custom*, *loca*, and *xyz* * element = "C" for every atom type * every = whatever it was set to via the :doc:`dump ` command diff --git a/src/dump.cpp b/src/dump.cpp index cdf2f57d9e..21803d79ab 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -1270,7 +1270,8 @@ void Dump::modify_params(int narg, char **arg) } else if (strcmp(arg[iarg],"balance") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); - balance_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + if (nprocs > 1) + balance_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"time") == 0) { From aeb25b5a3729ed2c5ecb43d8897ea04794a75fc0 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 15:28:40 -0700 Subject: [PATCH 64/75] Tweak docs --- doc/src/dump_modify.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index 52efa47d7c..5f77b79347 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -661,10 +661,11 @@ order. In a parallel run, the per-processor dump file pieces can have significant imbalance in number of lines of per-atom info. The *balance* -keyword determines whether lines of per-atom output for each processor -snapshot in a parallel run are balanced to be nearly the same. A balance -value of *no* means no balancing will be done, while *yes* means -balancing will be performed. For a serial run, this option is ignored. +keyword determines whether the number of lines in each processor +snapshot are balanced to be nearly the same. A balance value of *no* +means no balancing will be done, while *yes* means balancing will be +performed. For a serial run, this option is ignored since the output is +already balanced. Dump sorting must be enabled to use balancing. The dump *local* style cannot be sorted by atom ID, since there are typically multiple lines of output per atom. Some dump styles, such From a80f7ed11fdce58dd2b2b38597ef7f1a91531490 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 15:32:20 -0700 Subject: [PATCH 65/75] Add error message to header file --- src/dump.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/dump.h b/src/dump.h index 83edda3810..9054546a2d 100644 --- a/src/dump.h +++ b/src/dump.h @@ -176,6 +176,10 @@ E: Dump file MPI-IO output not allowed with % in filename This is because a % signifies one file per processor and MPI-IO creates one large file for all processors. +E: Cannot balance dump output without sorting enabled + +Self-explanatory. + E: Cannot dump sort when 'nfile' or 'fileper' keywords are set to non-default values Can only dump sort when the number of dump file pieces using % in filename equals the number of processors From cb9ba8c0a3fb7b4c3876cea64f80b9e3003abe34 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 16:05:27 -0700 Subject: [PATCH 66/75] Simplify code --- src/dump.cpp | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/dump.cpp b/src/dump.cpp index 21803d79ab..6bc7f401a6 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -937,9 +937,6 @@ void Dump::balance() proc_new_offsets[0] = 0; - bigint start = proc_new_offsets[me]; - bigint end = proc_new_offsets[me+1]; - // reset buf size to largest of any post-balance nme values // this insures proc 0 can receive everyone's info @@ -956,18 +953,18 @@ void Dump::balance() int nswap = 0; MPI_Request *request = new MPI_Request[nprocs]; - int procstart = start; + int procstart = 0; int iproc = me; int iproc_prev; - for (bigint i = start; i < end; i++) { + for (int i = 0; i < nme_balance; i++) { // find which proc this atom belongs to - while (i < proc_offsets[iproc]) iproc--; - while (i > proc_offsets[iproc+1]-1) iproc++; + while (proc_new_offsets[me] + i < proc_offsets[iproc]) iproc--; + while (proc_new_offsets[me] + i > proc_offsets[iproc+1]-1) iproc++; - if (i != start && (iproc != iproc_prev || i == end-1)) { + if (i != 0 && (iproc != iproc_prev || i == nme_balance - 1)) { // finished with proc @@ -980,7 +977,7 @@ void Dump::balance() // post receive for this proc if (iproc_prev != me) - MPI_Irecv(&buf_balance[(procstart - start)*size_one],procnrecv*size_one,MPI_DOUBLE, + MPI_Irecv(&buf_balance[procstart*size_one],procnrecv*size_one,MPI_DOUBLE, procrecv,0,world,&request[nswap++]); procstart = i; @@ -1018,7 +1015,7 @@ void Dump::balance() // sending to self, copy buffers - int offset_me = proc_offsets[me] - start; + int offset_me = proc_offsets[me] - proc_new_offsets[me]; memcpy(&buf_balance[(offset_me + procstart)*size_one],&buf[procstart*size_one],procnsend*size_one*sizeof(double)); } From 67af170929006eb73e9108e5fa27b0f22902d9d6 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 16:22:55 -0700 Subject: [PATCH 67/75] Sorting is not required for balancing --- doc/src/dump_modify.rst | 16 ++++++++-------- src/dump.cpp | 3 --- src/dump.h | 4 ---- 3 files changed, 8 insertions(+), 15 deletions(-) diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index 5f77b79347..d3a79ab073 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -659,14 +659,6 @@ atom ID. A sort value of N or -N means sort the output by the value in the Nth column of per-atom info in either ascending or descending order. -In a parallel run, the per-processor dump file pieces can have -significant imbalance in number of lines of per-atom info. The *balance* -keyword determines whether the number of lines in each processor -snapshot are balanced to be nearly the same. A balance value of *no* -means no balancing will be done, while *yes* means balancing will be -performed. For a serial run, this option is ignored since the output is -already balanced. Dump sorting must be enabled to use balancing. - The dump *local* style cannot be sorted by atom ID, since there are typically multiple lines of output per atom. Some dump styles, such as *dcd* and *xtc*, require sorting by atom ID to format the output @@ -676,6 +668,14 @@ keywords are set to non-default values (i.e. the number of dump file pieces is not equal to the number of procs), then sorting cannot be performed. +In a parallel run, the per-processor dump file pieces can have +significant imbalance in number of lines of per-atom info. The *balance* +keyword determines whether the number of lines in each processor +snapshot are balanced to be nearly the same. A balance value of *no* +means no balancing will be done, while *yes* means balancing will be +performed. For a serial run, this option is ignored since the output is +already balanced. + .. note:: Unless it is required by the dump style, sorting dump file diff --git a/src/dump.cpp b/src/dump.cpp index 6bc7f401a6..1043970d33 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -223,9 +223,6 @@ void Dump::init() ids = idsort = nullptr; index = proclist = nullptr; irregular = nullptr; - - if (balance_flag) - error->all(FLERR,"Cannot balance dump output without sorting enabled"); } if (sort_flag) { diff --git a/src/dump.h b/src/dump.h index 9054546a2d..83edda3810 100644 --- a/src/dump.h +++ b/src/dump.h @@ -176,10 +176,6 @@ E: Dump file MPI-IO output not allowed with % in filename This is because a % signifies one file per processor and MPI-IO creates one large file for all processors. -E: Cannot balance dump output without sorting enabled - -Self-explanatory. - E: Cannot dump sort when 'nfile' or 'fileper' keywords are set to non-default values Can only dump sort when the number of dump file pieces using % in filename equals the number of processors From c914fc6bafa2764a990bc6f023b6bb103e997f59 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 16:33:11 -0700 Subject: [PATCH 68/75] Small doc tweak --- doc/src/dump_modify.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index d3a79ab073..b25e5815a7 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -673,8 +673,8 @@ significant imbalance in number of lines of per-atom info. The *balance* keyword determines whether the number of lines in each processor snapshot are balanced to be nearly the same. A balance value of *no* means no balancing will be done, while *yes* means balancing will be -performed. For a serial run, this option is ignored since the output is -already balanced. +performed. This balancing preserves dump sorting order. For a serial +run, this option is ignored since the output is already balanced. .. note:: From 0b693e729df89ad575c32894b36876f9930f6f6b Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 24 Jan 2022 17:33:36 -0700 Subject: [PATCH 69/75] Whitespace --- doc/src/dump_modify.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index b25e5815a7..352f9c61bf 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -674,7 +674,7 @@ keyword determines whether the number of lines in each processor snapshot are balanced to be nearly the same. A balance value of *no* means no balancing will be done, while *yes* means balancing will be performed. This balancing preserves dump sorting order. For a serial -run, this option is ignored since the output is already balanced. +run, this option is ignored since the output is already balanced. .. note:: From 35dbabd4716765aa8886d01c7384a24a1393e3e9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 25 Jan 2022 07:39:37 -0500 Subject: [PATCH 70/75] alternate fix for tag caching issue in INTEL package --- src/INTEL/fix_intel.cpp | 2 +- src/INTEL/intel_buffers.cpp | 2 -- src/INTEL/npair_intel.cpp | 4 +++- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/INTEL/fix_intel.cpp b/src/INTEL/fix_intel.cpp index 6a848a366e..a201386724 100644 --- a/src/INTEL/fix_intel.cpp +++ b/src/INTEL/fix_intel.cpp @@ -447,7 +447,7 @@ void FixIntel::pair_init_check(const bool cdmessage) #endif int need_tag = 0; - if (atom->molecular != Atom::ATOMIC) need_tag = 1; + if (atom->molecular != Atom::ATOMIC || three_body_neighbor()) need_tag = 1; // Clear buffers used for pair style char kmode[80]; diff --git a/src/INTEL/intel_buffers.cpp b/src/INTEL/intel_buffers.cpp index 3470c18db8..cbbc609fc0 100644 --- a/src/INTEL/intel_buffers.cpp +++ b/src/INTEL/intel_buffers.cpp @@ -207,8 +207,6 @@ void IntelBuffers::free_nmax() template void IntelBuffers::_grow_nmax(const int offload_end) { - if (lmp->atom->molecular) _need_tag = 1; - else _need_tag = 0; #ifdef _LMP_INTEL_OFFLOAD free_nmax(); int size = lmp->atom->nmax; diff --git a/src/INTEL/npair_intel.cpp b/src/INTEL/npair_intel.cpp index 67cd2517e4..2fffb2353a 100644 --- a/src/INTEL/npair_intel.cpp +++ b/src/INTEL/npair_intel.cpp @@ -309,6 +309,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, flt_t * _noalias const tz = ncachez + toffs; int * _noalias const tj = ncachej + toffs; int * _noalias const tjtype = ncachejtype + toffs; + tagint * _noalias const ttag = ncachetag + toffs; flt_t * _noalias itx; flt_t * _noalias ity; @@ -365,6 +366,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, ty[u] = x[j].y; tz[u] = x[j].z; tjtype[u] = x[j].w; + if (THREE) ttag[u] = tag[j]; } if (FULL == 0 && TRI != 1) { @@ -512,7 +514,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, } if (THREE) { - const tagint jtag = tag[j]; + const tagint jtag = ttag[u]; int flist = 0; if (itag > jtag) { if (((itag+jtag) & 1) == 0) flist = 1; From 64dbce26dce54fe77839df0d720bea3f493fd348 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Tue, 25 Jan 2022 08:45:13 -0700 Subject: [PATCH 71/75] Cannot shrink buf --- src/dump.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/dump.cpp b/src/dump.cpp index 1043970d33..732f4a73ef 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -939,11 +939,14 @@ void Dump::balance() int nmax; MPI_Allreduce(&nme_balance,&nmax,1,MPI_INT,MPI_MAX,world); + if (nmax > maxbuf) { + maxbuf = nmax; + } // allocate a second buffer for balanced data double* buf_balance; - memory->create(buf_balance,nmax*size_one,"dump:buf_balance"); + memory->create(buf_balance,maxbuf*size_one,"dump:buf_balance"); // compute from which procs I am receiving atoms // post recvs first From 75a00f4d11a538d99509fd71326a51c1bb3dc94b Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Tue, 25 Jan 2022 08:52:53 -0700 Subject: [PATCH 72/75] minor cleanup --- src/dump.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/dump.cpp b/src/dump.cpp index 732f4a73ef..181a56ff46 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -936,12 +936,11 @@ void Dump::balance() // reset buf size to largest of any post-balance nme values // this insures proc 0 can receive everyone's info + // cannot shrink buf to nme_balance, must use previous maxbuf value int nmax; MPI_Allreduce(&nme_balance,&nmax,1,MPI_INT,MPI_MAX,world); - if (nmax > maxbuf) { - maxbuf = nmax; - } + if (nmax > maxbuf) maxbuf = nmax; // allocate a second buffer for balanced data From 120662c438be8bb5fa8da37543426c7458ac3945 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 26 Jan 2022 10:45:30 -0500 Subject: [PATCH 73/75] correct citation info --- src/MC/fix_bond_swap.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 5f7ffcbbe5..c0978f51cb 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -39,7 +39,7 @@ using namespace LAMMPS_NS; using namespace FixConst; static const char cite_fix_bond_swap[] = - "neighbor multi command:\n\n" + "fix bond/swap command:\n\n" "@Article{Auhl03,\n" " author = {R. Auhl, R. Everaers, G. S. Grest, K. Kremer, S. J. Plimpton},\n" " title = {Equilibration of long chain polymer melts in computer simulations},\n" From 8ca1004c03d527b6ba827436161ae45ddd45fb11 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 26 Jan 2022 10:45:46 -0500 Subject: [PATCH 74/75] cosmetic. align with clang-format --- src/RIGID/fix_rigid_small.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 14742155db..cbaef073da 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -220,7 +220,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"infile") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); - delete [] inpfile; + delete[] inpfile; inpfile = utils::strdup(arg[iarg+1]); restart_file = 1; reinitflag = 0; @@ -327,7 +327,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; else { allremap = 0; - delete [] id_dilate; + delete[] id_dilate; id_dilate = utils::strdup(arg[iarg+1]); int idilate = group->find(id_dilate); if (idilate == -1) @@ -354,7 +354,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"gravity") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); - delete [] id_gravity; + delete[] id_gravity; id_gravity = utils::strdup(arg[iarg+1]); iarg += 2; @@ -395,7 +395,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : double time1 = platform::walltime(); create_bodies(bodyID); - if (customflag) delete [] bodyID; + if (customflag) delete[] bodyID; if (comm->me == 0) utils::logmesg(lmp," create bodies CPU = {:.3f} seconds\n", @@ -501,9 +501,9 @@ FixRigidSmall::~FixRigidSmall() memory->destroy(dorient); delete random; - delete [] inpfile; - delete [] id_dilate; - delete [] id_gravity; + delete[] inpfile; + delete[] id_dilate; + delete[] id_gravity; memory->destroy(langextra); memory->destroy(mass_body); @@ -2578,8 +2578,8 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody) if (me == 0) fclose(fp); - delete [] buffer; - delete [] values; + delete[] buffer; + delete[] values; } /* ---------------------------------------------------------------------- From 71a373cadbf22214c3400d4312a72e59f9d07a2e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 26 Jan 2022 10:46:19 -0500 Subject: [PATCH 75/75] reshuffle Body struct members and add dummy for better alignment --- src/RIGID/fix_rigid_small.h | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index e289c179d9..27076ccac3 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -84,8 +84,9 @@ class FixRigidSmall : public Fix { double maxextent; // furthest distance from body owner to body atom struct Body { - double mass; // total mass of body int natoms; // total number of atoms in body + int ilocal; // index of owning atom + double mass; // total mass of body double xcm[3]; // COM position double xgc[3]; // geometric center position double vcm[3]; // COM velocity @@ -97,12 +98,12 @@ class FixRigidSmall : public Fix { double ey_space[3]; double ez_space[3]; double xgc_body[3]; // geometric center relative to xcm in body coords - double angmom[3]; // space-frame angular momentum of body - double omega[3]; // space-frame omega of body - double conjqm[4]; // conjugate quaternion momentum - imageint image; // image flags of xcm - int remapflag[4]; // PBC remap flags - int ilocal; // index of owning atom + double angmom[3]; // space-frame angular momentum of body + double omega[3]; // space-frame omega of body + double conjqm[4]; // conjugate quaternion momentum + int remapflag[4]; // PBC remap flags + imageint image; // image flags of xcm + imageint dummy; // dummy entry for better alignment }; Body *body; // list of rigid bodies, owned and ghost