From ba27ecf6102b17d1f92037e488effbcad279cd77 Mon Sep 17 00:00:00 2001 From: rohskopf Date: Mon, 3 Oct 2022 11:53:41 -0600 Subject: [PATCH 01/21] Add optional full flag to pair_zero --- src/pair_zero.cpp | 30 ++++++++++++++++++++++++++++-- src/pair_zero.h | 2 ++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/pair_zero.cpp b/src/pair_zero.cpp index 83505edd35..e253cd2aa5 100644 --- a/src/pair_zero.cpp +++ b/src/pair_zero.cpp @@ -21,6 +21,8 @@ #include "comm.h" #include "error.h" #include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" #include @@ -34,6 +36,7 @@ PairZero::PairZero(LAMMPS *lmp) : Pair(lmp) writedata = 1; single_enable = 1; respa_enable = 1; + fullneighflag = 0; } /* ---------------------------------------------------------------------- */ @@ -85,12 +88,23 @@ void PairZero::allocate() void PairZero::settings(int narg, char **arg) { - if ((narg != 1) && (narg != 2)) error->all(FLERR, "Illegal pair_style command"); + if ((narg != 1) && (narg != 2) && (narg != 3)) error->all(FLERR, "Illegal pair_style command"); cut_global = utils::numeric(FLERR, arg[0], false, lmp); - if (narg == 2) { + + if (narg == 2){ if (strcmp("nocoeff", arg[1]) == 0) coeffflag = 0; + else if(strcmp("full", arg[1]) == 0) + fullneighflag = 1; + else + error->all(FLERR, "Illegal pair_style command"); + } + else if (narg == 3) { + if (strcmp("nocoeff", arg[1]) == 0 || strcmp("nocoeff", arg[2]) == 0) + coeffflag = 0; + else if(strcmp("full", arg[1]) == 0 || strcmp("full", arg[2]) == 0) + fullneighflag = 1; else error->all(FLERR, "Illegal pair_style command"); } @@ -134,6 +148,18 @@ void PairZero::coeff(int narg, char **arg) if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairZero::init_style() +{ + if (fullneighflag == 0) + neighbor->add_request(this, NeighConst::REQ_DEFAULT); + else + neighbor->add_request(this, NeighConst::REQ_FULL); +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/pair_zero.h b/src/pair_zero.h index 03b6d096a1..10102f1cad 100644 --- a/src/pair_zero.h +++ b/src/pair_zero.h @@ -42,6 +42,7 @@ class PairZero : public Pair { void compute_outer(int, int) override; void settings(int, char **) override; void coeff(int, char **) override; + void init_style() override; double init_one(int, int) override; void write_restart(FILE *) override; void read_restart(FILE *) override; @@ -55,6 +56,7 @@ class PairZero : public Pair { double cut_global; double **cut; int coeffflag; + int fullneighflag; // 0 for half list, 1 for full list virtual void allocate(); }; From c82fdd4898719dadca4ffc72746e1a5acaeefa3b Mon Sep 17 00:00:00 2001 From: rohskopf Date: Mon, 3 Oct 2022 11:57:47 -0600 Subject: [PATCH 02/21] Add docs --- doc/src/pair_zero.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/src/pair_zero.rst b/doc/src/pair_zero.rst index 062b5204f8..0e6bac643a 100644 --- a/doc/src/pair_zero.rst +++ b/doc/src/pair_zero.rst @@ -8,11 +8,12 @@ Syntax .. code-block:: LAMMPS - pair_style zero cutoff [nocoeff] + pair_style zero cutoff [nocoeff] [full] * zero = style name of this pair style * cutoff = global cutoff (distance units) * nocoeff = ignore all pair_coeff parameters (optional) +* full = build full neighbor list (optional) Examples """""""" @@ -45,6 +46,9 @@ section for any pair style. Similarly, any pair_coeff commands will only be checked for the atom type numbers and the rest ignored. In this case, only the global cutoff will be used. +The optional *full* flag builds a full neighbor list instead of the default +half neighbor list. + 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 From 981f75d92fa685b25d5f2828dc3c3635910fa8c8 Mon Sep 17 00:00:00 2001 From: rohskopf Date: Mon, 3 Oct 2022 12:40:05 -0600 Subject: [PATCH 03/21] Example in docs --- doc/src/Python_neighbor.rst | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/doc/src/Python_neighbor.rst b/doc/src/Python_neighbor.rst index 00c4cc8996..8201955b6c 100644 --- a/doc/src/Python_neighbor.rst +++ b/doc/src/Python_neighbor.rst @@ -38,6 +38,40 @@ using the NumPy access method. for n in np.nditer(nlist): print(" atom {} with ID {}".format(n,tags[n])) +Another example for extracting a full neighbor list without evaluating a +potential is shown below. + +.. code-block:: python + + from lammps import lammps + import numpy as np + + lmp = lammps() + lmp.commands_string(""" + newton off + region box block -2 2 -2 2 -2 2 + lattice fcc 1.0 + create_box 1 box + create_atoms 1 box + mass 1 1.0 + pair_style zero 1.0 full + pair_coeff * * + run 0 post no""") + + # look up the neighbor list + nlidx = lmp.find_pair_neighlist('zero') + nl = lmp.numpy.get_neighlist(nlidx) + tags = lmp.extract_atom('id') + print("full neighbor list with {} entries".format(nl.size)) + # print neighbor list contents + for i in range(0,nl.size): + idx, nlist = nl.get(i) + print("\natom {} with ID {} has {} neighbors:".format(idx,tags[idx],nlist.size)) + if nlist.size > 0: + for n in np.nditer(nlist): + pass + print(" atom {} with ID {}".format(n,tags[n])) + **Methods:** * :py:meth:`lammps.get_neighlist() `: Get neighbor list for given index From 82bf23401f4747338eacc3d5c2d94ece012d8ab3 Mon Sep 17 00:00:00 2001 From: rohskopf Date: Mon, 3 Oct 2022 12:45:18 -0600 Subject: [PATCH 04/21] Fix whitespace --- doc/src/Python_neighbor.rst | 2 +- src/pair_zero.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/Python_neighbor.rst b/doc/src/Python_neighbor.rst index 8201955b6c..755b83da81 100644 --- a/doc/src/Python_neighbor.rst +++ b/doc/src/Python_neighbor.rst @@ -38,7 +38,7 @@ using the NumPy access method. for n in np.nditer(nlist): print(" atom {} with ID {}".format(n,tags[n])) -Another example for extracting a full neighbor list without evaluating a +Another example for extracting a full neighbor list without evaluating a potential is shown below. .. code-block:: python diff --git a/src/pair_zero.cpp b/src/pair_zero.cpp index e253cd2aa5..e4d98ad66e 100644 --- a/src/pair_zero.cpp +++ b/src/pair_zero.cpp @@ -91,7 +91,7 @@ void PairZero::settings(int narg, char **arg) if ((narg != 1) && (narg != 2) && (narg != 3)) error->all(FLERR, "Illegal pair_style command"); cut_global = utils::numeric(FLERR, arg[0], false, lmp); - + if (narg == 2){ if (strcmp("nocoeff", arg[1]) == 0) coeffflag = 0; From 21aded5e4ee6139e9ddefecf8f729f210251cc5b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Oct 2022 15:14:54 -0400 Subject: [PATCH 05/21] silence compiler warning --- src/fix_ave_time.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index 3eeed90e2f..bd8ed1468d 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -895,7 +895,7 @@ void FixAveTime::invoke_vector(bigint ntimestep) int FixAveTime::column_length(int dynamic) { - int m,length,lengthone; + int length,lengthone; // determine nrows for static values From 840e3398b8b99aadb5cbe104c1c7a6ce2b323956 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Oct 2022 15:36:19 -0400 Subject: [PATCH 06/21] silence warning with CMake 3.24 and later --- cmake/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 28c3d6c027..8267a9bf91 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -12,6 +12,11 @@ endif() if(POLICY CMP0109) cmake_policy(SET CMP0109 OLD) endif() +# set policy to silence warnings about timestamps of downloaded files. review occasionally if it may be set to NEW +if(POLICY CMP0135) + cmake_policy(SET CMP0135 OLD) +endif() + ######################################## project(lammps CXX) From 2d3e4fdd9a601868cb4e77db4073ac374b9392b1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Oct 2022 15:24:05 -0400 Subject: [PATCH 07/21] make MDI package cmake config compatible with multi-config builders --- cmake/Modules/Packages/MDI.cmake | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/cmake/Modules/Packages/MDI.cmake b/cmake/Modules/Packages/MDI.cmake index edae4ffcbd..82b036a583 100644 --- a/cmake/Modules/Packages/MDI.cmake +++ b/cmake/Modules/Packages/MDI.cmake @@ -70,23 +70,22 @@ if(DOWNLOAD_MDI) -Dplugins=ON -Dpython_plugins=${MDI_USE_PYTHON_PLUGINS} UPDATE_COMMAND "" - INSTALL_COMMAND "" - BUILD_BYPRODUCTS "/MDI_Library/libmdi.a" + BUILD_BYPRODUCTS "/${CMAKE_INSTALL_LIBDIR}/mdi/${CMAKE_STATIC_LIBRARY_PREFIX}mdi${CMAKE_STATIC_LIBRARY_SUFFIX}" ) # where is the compiled library? - ExternalProject_get_property(mdi_build BINARY_DIR) - set(MDI_BINARY_DIR "${BINARY_DIR}/MDI_Library") + ExternalProject_get_property(mdi_build INSTALL_DIR) # workaround for older CMake versions - file(MAKE_DIRECTORY ${MDI_BINARY_DIR}) + file(MAKE_DIRECTORY ${INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/mdi) + file(MAKE_DIRECTORY ${INSTALL_DIR}/include/mdi) # create imported target for the MDI library add_library(LAMMPS::MDI UNKNOWN IMPORTED) add_dependencies(LAMMPS::MDI mdi_build) set_target_properties(LAMMPS::MDI PROPERTIES - IMPORTED_LOCATION "${MDI_BINARY_DIR}/libmdi.a" - INTERFACE_INCLUDE_DIRECTORIES ${MDI_BINARY_DIR} - ) + IMPORTED_LOCATION "${INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/mdi/${CMAKE_STATIC_LIBRARY_PREFIX}mdi${CMAKE_STATIC_LIBRARY_SUFFIX}" + INTERFACE_INCLUDE_DIRECTORIES ${INSTALL_DIR}/include/mdi + ) set(MDI_DEP_LIBS "") # if compiling with python plugins we need From bb6a1180067e6de2639120ff5742ff9dfb492f96 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Oct 2022 16:31:23 -0400 Subject: [PATCH 08/21] refactor for simpler flow of control. apply clang-format. --- src/pair_zero.cpp | 33 ++++++++++++++++----------------- src/pair_zero.h | 4 +--- 2 files changed, 17 insertions(+), 20 deletions(-) diff --git a/src/pair_zero.cpp b/src/pair_zero.cpp index e4d98ad66e..b554f78b7f 100644 --- a/src/pair_zero.cpp +++ b/src/pair_zero.cpp @@ -88,25 +88,24 @@ void PairZero::allocate() void PairZero::settings(int narg, char **arg) { - if ((narg != 1) && (narg != 2) && (narg != 3)) error->all(FLERR, "Illegal pair_style command"); + if (narg < 1) utils::missing_cmd_args(FLERR, "pair_style zero", error); cut_global = utils::numeric(FLERR, arg[0], false, lmp); - if (narg == 2){ - if (strcmp("nocoeff", arg[1]) == 0) + // reset to defaults + coeffflag = 1; + fullneighflag = 0; + + int iarg = 1; + while (iarg < narg) { + if (strcmp("nocoeff", arg[iarg]) == 0) { coeffflag = 0; - else if(strcmp("full", arg[1]) == 0) + ++iarg; + } else if (strcmp("full", arg[iarg]) == 0) { fullneighflag = 1; - else - error->all(FLERR, "Illegal pair_style command"); - } - else if (narg == 3) { - if (strcmp("nocoeff", arg[1]) == 0 || strcmp("nocoeff", arg[2]) == 0) - coeffflag = 0; - else if(strcmp("full", arg[1]) == 0 || strcmp("full", arg[2]) == 0) - fullneighflag = 1; - else - error->all(FLERR, "Illegal pair_style command"); + ++iarg; + } else + error->all(FLERR, "Unknown pair style zero option {}", arg[iarg]); } // reset cutoffs that have been explicitly set @@ -154,10 +153,10 @@ void PairZero::coeff(int narg, char **arg) void PairZero::init_style() { - if (fullneighflag == 0) - neighbor->add_request(this, NeighConst::REQ_DEFAULT); - else + if (fullneighflag) neighbor->add_request(this, NeighConst::REQ_FULL); + else + neighbor->add_request(this, NeighConst::REQ_DEFAULT); } /* ---------------------------------------------------------------------- diff --git a/src/pair_zero.h b/src/pair_zero.h index 10102f1cad..06d7b2ab72 100644 --- a/src/pair_zero.h +++ b/src/pair_zero.h @@ -56,12 +56,10 @@ class PairZero : public Pair { double cut_global; double **cut; int coeffflag; - int fullneighflag; // 0 for half list, 1 for full list + int fullneighflag; // 0 for half list, 1 for full list virtual void allocate(); }; - } // namespace LAMMPS_NS - #endif #endif From ef3b01f3402cb92ec20902c7c469fdbde2e8b2e9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Oct 2022 18:14:08 -0400 Subject: [PATCH 09/21] add neighbor list access tests. now we can always test for a full list. --- unittest/python/python-commands.py | 69 ++++++++++++++++++++++++++++++ unittest/python/python-numpy.py | 68 +++++++++++++++++++++++++++++ 2 files changed, 137 insertions(+) diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py index cd77fa21a1..13be27f067 100644 --- a/unittest/python/python-commands.py +++ b/unittest/python/python-commands.py @@ -209,6 +209,75 @@ create_atoms 1 single & self.assertEqual(idx,i) self.assertEqual(num,nlocal-1) + def testNeighborListZeroHalf(self): + self.lmp.commands_string(""" + boundary f f f + units real + region box block -5 5 -5 5 -5 5 + create_box 1 box + mass 1 1.0 + pair_style zero 4.0 + pair_coeff 1 1 + """) + x = [ 0.0, 0.0, 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.1, + 0.0, 0.0, 1.0 ] + tags = [1, 2, 3, 4, 5, 6, 7] + types = [1, 1, 1, 1, 1, 1, 1] + + self.assertEqual(self.lmp.create_atoms(7, id=tags, type=types, x=x), 7) + nlocal = self.lmp.extract_global("nlocal") + self.assertEqual(nlocal, 7) + + self.lmp.command("run 0 post no") + + self.assertEqual(self.lmp.find_pair_neighlist("zero"),0) + nlist = self.lmp.get_neighlist(0) + self.assertEqual(nlist.size, 7) + for i in range(0,nlist.size): + idx, num, neighs = nlist.get(i) + self.assertEqual(idx,i) + self.assertEqual(num,nlocal-1-i) + + # look up neighbor list by atom index + num, neighs = nlist.find(2) + self.assertEqual(num,4) + self.assertIsNotNone(neighs,None) + # this one will fail + num, neighs = nlist.find(10) + self.assertEqual(num,-1) + self.assertIsNone(neighs,None) + + def testNeighborListZeroFull(self): + self.lmp.commands_string(""" + boundary f f f + units metal + region box block -5 5 -5 5 -5 5 + create_box 1 box + mass 1 1.0 + pair_style zero 4.0 full + pair_coeff * * + """) + x = [ 0.0, 0.0, 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.1, + 0.0, 0.0, 1.0 ] + tags = [1, 2, 3, 4, 5, 6, 7] + types = [1, 1, 1, 1, 1, 1, 1] + + self.assertEqual(self.lmp.create_atoms(7, id=tags, type=types, x=x), 7) + nlocal = self.lmp.extract_global("nlocal") + self.assertEqual(nlocal, 7) + + self.lmp.command("run 0 post no") + + self.assertEqual(self.lmp.find_pair_neighlist("zero"),0) + nlist = self.lmp.get_neighlist(0) + self.assertEqual(nlist.size, 7) + for i in range(0,nlist.size): + idx, num, neighs = nlist.get(i) + self.assertEqual(idx,i) + self.assertEqual(num,nlocal-1) + @unittest.skipIf(not has_manybody,"Hybrid neighbor list test for manybody potential") def testNeighborListHybrid(self): self.lmp.commands_string(""" diff --git a/unittest/python/python-numpy.py b/unittest/python/python-numpy.py index 7cbae8e48d..3306aba1cd 100644 --- a/unittest/python/python-numpy.py +++ b/unittest/python/python-numpy.py @@ -412,6 +412,74 @@ class PythonNumpy(unittest.TestCase): idx, neighs = nlist.get(i) self.assertEqual(neighs.size,nlocal-1-i) + def testNeighborListZeroHalf(self): + self.lmp.commands_string(""" + boundary f f f + units real + region box block -5 5 -5 5 -5 5 + create_box 1 box + mass 1 1.0 + pair_style zero 4.0 + pair_coeff 1 1 + """) + x = [ 0.0, 0.0, 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.1, + 0.0, 0.0, 1.0 ] + tags = [1, 2, 3, 4, 5, 6, 7] + types = [1, 1, 1, 1, 1, 1, 1] + + self.assertEqual(self.lmp.create_atoms(7, id=tags, type=types, x=x), 7) + nlocal = self.lmp.extract_global("nlocal") + self.assertEqual(nlocal, 7) + + self.lmp.command("run 0 post no") + + self.assertEqual(self.lmp.find_pair_neighlist("zero"),0) + nlist = self.lmp.numpy.get_neighlist(0) + self.assertEqual(nlist.size, 7) + for i in range(0,nlist.size): + idx, neighs = nlist.get(i) + self.assertEqual(idx,i) + self.assertEqual(neighs.size,nlocal-1-i) + + # look up neighbor list by atom index + neighs = nlist.find(2) + self.assertEqual(neighs.size,4) + self.assertIsNotNone(neighs,None) + # this one will fail + neighs = nlist.find(10) + self.assertIsNone(neighs,None) + + def testNeighborListZeroFull(self): + self.lmp.commands_string(""" + boundary f f f + units metal + region box block -5 5 -5 5 -5 5 + create_box 1 box + mass 1 1.0 + pair_style zero 4.0 full + pair_coeff * * + """) + x = [ 0.0, 0.0, 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, -1.1, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.1, + 0.0, 0.0, 1.0 ] + tags = [1, 2, 3, 4, 5, 6, 7] + types = [1, 1, 1, 1, 1, 1, 1] + + self.assertEqual(self.lmp.create_atoms(7, id=tags, type=types, x=x), 7) + nlocal = self.lmp.extract_global("nlocal") + self.assertEqual(nlocal, 7) + + self.lmp.command("run 0 post no") + + self.assertEqual(self.lmp.find_pair_neighlist("zero"),0) + nlist = self.lmp.numpy.get_neighlist(0) + self.assertEqual(nlist.size, 7) + for i in range(0,nlist.size): + idx, neighs = nlist.get(i) + self.assertEqual(idx,i) + self.assertEqual(neighs.size,nlocal-1) + def testNeighborListCompute(self): self.lmp.commands_string(""" boundary f f f From ed8e317ef62646c01c0e4b90253f8f85810f9005 Mon Sep 17 00:00:00 2001 From: Terry Suun Date: Tue, 4 Oct 2022 10:40:32 +0800 Subject: [PATCH 10/21] Added quartic function for explicit pair interaction in pair_style list. --- doc/src/pair_list.rst | 16 ++++++++++++++++ src/MISC/pair_list.cpp | 42 +++++++++++++++++++++++++++++++++++------- src/MISC/pair_list.h | 5 +++++ 3 files changed, 56 insertions(+), 7 deletions(-) diff --git a/doc/src/pair_list.rst b/doc/src/pair_list.rst index 321bf61be9..322da79369 100644 --- a/doc/src/pair_list.rst +++ b/doc/src/pair_list.rst @@ -69,6 +69,7 @@ Here is an example file: 15 259 lj126 1.0 1.0 50.0 15 603 morse 10.0 1.2 2.0 10.0 # and another comment 18 470 harmonic 50.0 1.2 5.0 + 19 332 quartic 5.0 -1.2 1.2 5.0 The style *lj126* computes pairwise interactions with the formula @@ -106,6 +107,21 @@ and the coefficients: Note that the usual 1/2 factor is included in :math:`K`. +The style *quartic* computes pairwise interactions with the formula + +.. math:: + + E = K (r - r_c)^2 (r - r_c -b_1) (r - r_c - b_2) \qquad r < r_c + +and the coefficients: + +* :math:`K` (energy units) +* :math:`b_1` (distance units) +* :math:`b_2` (distance units) +* :math:`r_c` (distance units) + +Note that the cutoff :math:`r_c` should always be specified to ensure zero energy and smooth force at cutoff. + ---------- Mixing, shift, table, tail correction, restart, rRESPA info diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index f6ac2e3190..1eff74564d 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -32,13 +32,14 @@ using namespace LAMMPS_NS; -enum { NONE = 0, HARM, MORSE, LJ126 }; +enum { NONE = 0, HARM, MORSE, LJ126, QUARTIC }; static std::map stylename = { { "none", NONE }, { "harmonic", HARM }, { "morse", MORSE }, - { "lj126", LJ126 } + { "lj126", LJ126 }, + { "quartic", QUARTIC } }; // fast power function for integer exponent > 0 @@ -157,8 +158,21 @@ void PairList::compute(int eflag, int vflag) if (eflag_either) epair = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6) - par.offset; + + } else if (par.style == QUARTIC) { + + const double r = sqrt(rsq); + double dr = r - sqrt(par.cutsq); + double ra = dr - par.param.quartic.b1; + double rb = dr - par.param.quartic.b2; + double r2 = dr * dr; + fpair = -par.param.quartic.k / r * (r2 * (ra + rb) + 2.0 * dr * ra * rb); + + if (eflag_either) + epair = par.param.quartic.k * r2 * ra * rb; } + if (newton_pair || i < nlocal) { f[i].x += dx*fpair; f[i].y += dy*fpair; @@ -221,10 +235,10 @@ void PairList::settings(int narg, char **arg) // read and parse potential file only on MPI rank 0. if (comm->me == 0) { - int nharm, nmorse, nlj126, nskipped; + int nharm, nmorse, nlj126, nquartic, nskipped; FILE *fp = utils::open_potential(arg[0],lmp,nullptr); TextFileReader reader(fp,"pair list coeffs"); - npairs = nharm = nmorse = nlj126 = nskipped = 0; + npairs = nharm = nmorse = nlj126 = nquartic = nskipped = 0; try { char *line; @@ -258,6 +272,14 @@ void PairList::settings(int narg, char **arg) ++nlj126; break; + case QUARTIC: + oneparam.param.quartic.k = values.next_double(); + oneparam.param.quartic.b1 = values.next_double(); + oneparam.param.quartic.b2 = values.next_double(); + oneparam.param.quartic.rc = values.next_double(); + ++nquartic; + break; + case NONE: // fallthrough error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}", utils::trim(line)); @@ -265,7 +287,7 @@ void PairList::settings(int narg, char **arg) break; } if (values.has_next()) - oneparam.cutsq = values.next_double(); + oneparam.cutsq = mypow(values.next_double(), 2); else oneparam.cutsq = cut_global*cut_global; @@ -274,8 +296,8 @@ void PairList::settings(int narg, char **arg) } 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); + utils::logmesg(lmp, "Read {} ({}/{}/{}/{}) interacting pair lines from {}. " + "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, nquartic, arg[0], nskipped); memory->create(params,npairs,"pair_list:params"); memcpy(params, myparams.data(),npairs*sizeof(list_param)); @@ -339,6 +361,12 @@ void PairList::init_style() const double r6inv = par.cutsq*par.cutsq*par.cutsq; const double sig6 = mypow(par.param.lj126.sigma,6); par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); + + } else if (par.style == QUARTIC) { + // the offset is always 0 at rc + par.offset = 0.0; + // correct cutsq + par.cutsq = mypow(par.param.quartic.rc, 2); } } } diff --git a/src/MISC/pair_list.h b/src/MISC/pair_list.h index 2828a08547..2d00788811 100644 --- a/src/MISC/pair_list.h +++ b/src/MISC/pair_list.h @@ -49,11 +49,16 @@ class PairList : public Pair { struct lj126_p { double epsilon, sigma; }; + struct quartic_p { + double k, b1, b2, rc; + }; + union param_u { harm_p harm; morse_p morse; lj126_p lj126; + quartic_p quartic; }; struct list_param { From bca8e6b85a06bc62b29b5712bd99cc1189db8898 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 00:14:06 -0400 Subject: [PATCH 11/21] store squared cutoff as documented --- src/MISC/pair_list.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index f6ac2e3190..be4e2cf92b 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -22,6 +22,7 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "math_special.h" #include "memory.h" #include "text_file_reader.h" @@ -31,6 +32,7 @@ #include using namespace LAMMPS_NS; +using MathSpecial::square; enum { NONE = 0, HARM, MORSE, LJ126 }; @@ -265,7 +267,7 @@ void PairList::settings(int narg, char **arg) break; } if (values.has_next()) - oneparam.cutsq = values.next_double(); + oneparam.cutsq = square(values.next_double()); else oneparam.cutsq = cut_global*cut_global; From 5508c7f25e2233960f3dbf139cc20129127e53c8 Mon Sep 17 00:00:00 2001 From: Terry Suun Date: Tue, 4 Oct 2022 13:27:32 +0800 Subject: [PATCH 12/21] Bug fix for pair_style list quartic. --- doc/src/pair_list.rst | 2 +- src/MISC/pair_list.cpp | 3 --- src/MISC/pair_list.h | 2 +- 3 files changed, 2 insertions(+), 5 deletions(-) diff --git a/doc/src/pair_list.rst b/doc/src/pair_list.rst index 322da79369..76b31623ad 100644 --- a/doc/src/pair_list.rst +++ b/doc/src/pair_list.rst @@ -120,7 +120,7 @@ and the coefficients: * :math:`b_2` (distance units) * :math:`r_c` (distance units) -Note that the cutoff :math:`r_c` should always be specified to ensure zero energy and smooth force at cutoff. +Note that the global cutoff specified by *pair_style list* command is ignored, unless :math:`r_c` is not specified. In this case, :math:`r_c` equals the sqare root of the globle cutoff. ---------- diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index 1eff74564d..ad3ef44717 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -276,7 +276,6 @@ void PairList::settings(int narg, char **arg) oneparam.param.quartic.k = values.next_double(); oneparam.param.quartic.b1 = values.next_double(); oneparam.param.quartic.b2 = values.next_double(); - oneparam.param.quartic.rc = values.next_double(); ++nquartic; break; @@ -365,8 +364,6 @@ void PairList::init_style() } else if (par.style == QUARTIC) { // the offset is always 0 at rc par.offset = 0.0; - // correct cutsq - par.cutsq = mypow(par.param.quartic.rc, 2); } } } diff --git a/src/MISC/pair_list.h b/src/MISC/pair_list.h index 2d00788811..dab8ba3fc9 100644 --- a/src/MISC/pair_list.h +++ b/src/MISC/pair_list.h @@ -50,7 +50,7 @@ class PairList : public Pair { double epsilon, sigma; }; struct quartic_p { - double k, b1, b2, rc; + double k, b1, b2; }; From a3c676015beb19aa100a1a0750c922dbd2c8c9bd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 04:13:42 -0400 Subject: [PATCH 13/21] store squared cutoff as documented --- src/MISC/pair_list.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index ad3ef44717..588af2b01a 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -22,6 +22,7 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "math_special.h" #include "memory.h" #include "text_file_reader.h" @@ -31,6 +32,7 @@ #include using namespace LAMMPS_NS; +using MathSpecial::square; enum { NONE = 0, HARM, MORSE, LJ126, QUARTIC }; @@ -286,7 +288,7 @@ void PairList::settings(int narg, char **arg) break; } if (values.has_next()) - oneparam.cutsq = mypow(values.next_double(), 2); + oneparam.cutsq = square(values.next_double()); else oneparam.cutsq = cut_global*cut_global; From 6a97ca246842096992d2a9cd6b0b0278be62309b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 05:02:46 -0400 Subject: [PATCH 14/21] better error checking and reporting --- src/MISC/pair_list.cpp | 82 +++++++++++++++++++++++++++--------------- 1 file changed, 54 insertions(+), 28 deletions(-) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index 588af2b01a..141afba4e6 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -89,7 +89,15 @@ void PairList::compute(int eflag, int vflag) { ev_init(eflag,vflag); + // get maximum allowed tag. + + bigint maxtag_one, maxtag; + maxtag_one = maxtag = 0; const int nlocal = atom->nlocal; + const tagint * _noalias const tag = atom->tag; + for (int i = 0; i < nlocal; ++i) maxtag_one = MAX(maxtag_one, tag[i]); + MPI_Allreduce(&maxtag_one, &maxtag, 1, MPI_LMP_TAGINT, MPI_MAX, world); + const int newton_pair = force->newton_pair; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; dbl3_t * _noalias const f = (dbl3_t *) atom->f[0]; // NOLINT @@ -98,8 +106,21 @@ void PairList::compute(int eflag, int vflag) int i,j; int pc = 0; - for (int n=0; n < npairs; ++n) { + for (int n = 0; n < npairs; ++n) { const list_param &par = params[n]; + + // can only use valid tags or else atom->map() below will segfault. + if ((par.id1 < 1) || (par.id1 > maxtag)) { + if (check_flag) + error->all(FLERR, "Invalid pair list atom ID {}", par.id1); + else continue; + } + if ((par.id2 < 1) || (par.id2 > maxtag)) { + if (check_flag) + error->all(FLERR, "Invalid pair list atom ID {}", par.id2); + else continue; + } + i = atom->map(par.id1); j = atom->map(par.id2); @@ -144,8 +165,7 @@ void PairList::compute(int eflag, int vflag) const double r = sqrt(rsq); 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; + fpair = 2.0 * par.param.morse.d0 * par.param.morse.alpha * (dexp*dexp - dexp) / r; if (eflag_either) epair = par.param.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; @@ -154,27 +174,24 @@ void PairList::compute(int eflag, int vflag) const double r6inv = r2inv*r2inv*r2inv; const double sig6 = mypow(par.param.lj126.sigma,6); - fpair = 24.0*par.param.lj126.epsilon*r6inv - * (2.0*sig6*sig6*r6inv - sig6) * r2inv; + fpair = 24.0 * par.param.lj126.epsilon * r6inv * (2.0 * sig6 * sig6 * r6inv - sig6) * r2inv; if (eflag_either) - epair = 4.0*par.param.lj126.epsilon*r6inv - * (sig6*sig6*r6inv - sig6) - par.offset; - + epair = 4.0 * par.param.lj126.epsilon * r6inv * (sig6 * sig6 * r6inv - sig6) - par.offset; + } else if (par.style == QUARTIC) { const double r = sqrt(rsq); double dr = r - sqrt(par.cutsq); - double ra = dr - par.param.quartic.b1; - double rb = dr - par.param.quartic.b2; - double r2 = dr * dr; + double ra = dr - par.param.quartic.b1; + double rb = dr - par.param.quartic.b2; + double r2 = dr * dr; fpair = -par.param.quartic.k / r * (r2 * (ra + rb) + 2.0 * dr * ra * rb); if (eflag_either) epair = par.param.quartic.k * r2 * ra * rb; } - if (newton_pair || i < nlocal) { f[i].x += dx*fpair; f[i].y += dy*fpair; @@ -207,14 +224,14 @@ void PairList::compute(int eflag, int vflag) void PairList::allocate() { allocated = 1; - int n = atom->ntypes; + int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) + memory->create(setflag,np1,np1,"pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq,np1,np1,"pair:cutsq"); } /* ---------------------------------------------------------------------- @@ -223,13 +240,18 @@ void PairList::allocate() void PairList::settings(int narg, char **arg) { - if (narg < 2) - error->all(FLERR,"Illegal pair_style command"); + if (narg < 2) utils::missing_cmd_args(FLERR, "pair_style list", error); cut_global = utils::numeric(FLERR,arg[1],false,lmp); - if (narg > 2) { - if (strcmp(arg[2],"nocheck") == 0) check_flag = 0; - if (strcmp(arg[2],"check") == 0) check_flag = 1; + int iarg = 2; + while (iarg < narg) { + if (strcmp(arg[iarg],"nocheck") == 0) { + check_flag = 0; + ++iarg; + } else if (strcmp(arg[2],"check") == 0) { + check_flag = 1; + ++iarg; + } else error->all(FLERR, "Unknown pair style list keyword: {}", arg[iarg]); } std::vector mystyles; @@ -239,11 +261,13 @@ void PairList::settings(int narg, char **arg) if (comm->me == 0) { int nharm, nmorse, nlj126, nquartic, nskipped; FILE *fp = utils::open_potential(arg[0],lmp,nullptr); - TextFileReader reader(fp,"pair list coeffs"); + if (!fp) + error->one(FLERR, "Error opening pair list coeffs file {}: {}", arg[0], utils::getsyserror()); + TextFileReader reader(fp, "pair list coeffs"); npairs = nharm = nmorse = nlj126 = nquartic = nskipped = 0; + char *line; try { - char *line; while ((line = reader.next_line())) { ValueTokenizer values(line); list_param oneparam; @@ -278,6 +302,8 @@ void PairList::settings(int narg, char **arg) oneparam.param.quartic.k = values.next_double(); oneparam.param.quartic.b1 = values.next_double(); oneparam.param.quartic.b2 = values.next_double(); + if (!values.has_next()) + throw FileReaderException("Must specify individual cutoff for quartic interaction"); ++nquartic; break; @@ -295,7 +321,7 @@ void PairList::settings(int narg, char **arg) myparams.push_back(oneparam); } } catch (std::exception &e) { - error->one(FLERR,"Error reading pair list coeffs file: {}", e.what()); + error->one(FLERR,"Error reading pair list coeffs file: {}\n{}", e.what(),line); } utils::logmesg(lmp, "Read {} ({}/{}/{}/{}) interacting pair lines from {}. " "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, nquartic, arg[0], nskipped); @@ -310,12 +336,12 @@ void PairList::settings(int narg, char **arg) } /* ---------------------------------------------------------------------- - there are no coeffs to be set, but we need to update setflag and pretend + there are no coeffs to be set, but we need to update setflag and pretend there are ------------------------------------------------------------------------- */ void PairList::coeff(int narg, char **arg) { - if (narg < 2) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 2) utils::missing_cmd_args(FLERR,"pair_coeff list", error); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -362,7 +388,7 @@ void PairList::init_style() const double r6inv = par.cutsq*par.cutsq*par.cutsq; const double sig6 = mypow(par.param.lj126.sigma,6); par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); - + } else if (par.style == QUARTIC) { // the offset is always 0 at rc par.offset = 0.0; From 74d5893dd937361a3b38e3753e12786f4c91d256 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 05:06:31 -0400 Subject: [PATCH 15/21] enable and apply clang-format --- src/MISC/pair_list.cpp | 221 +++++++++++++++++++++-------------------- 1 file changed, 112 insertions(+), 109 deletions(-) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index 141afba4e6..e826ff0693 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -26,8 +25,8 @@ #include "memory.h" #include "text_file_reader.h" -#include #include +#include #include #include @@ -36,27 +35,32 @@ using MathSpecial::square; enum { NONE = 0, HARM, MORSE, LJ126, QUARTIC }; +// clang-format off static std::map stylename = { - { "none", NONE }, - { "harmonic", HARM }, - { "morse", MORSE }, - { "lj126", LJ126 }, - { "quartic", QUARTIC } + {"none", NONE}, + {"harmonic", HARM}, + {"morse", MORSE}, + {"lj126", LJ126}, + {"quartic", QUARTIC} }; +// clang-format on // fast power function for integer exponent > 0 -static double mypow(double x, int n) { +static double mypow(double x, int n) +{ double yy; if (x == 0.0) return 0.0; - for (yy = 1.0; n != 0; n >>= 1, x *=x) + for (yy = 1.0; n != 0; n >>= 1, x *= x) if (n & 1) yy *= x; return yy; } -typedef struct { double x,y,z; } dbl3_t; +typedef struct { + double x, y, z; +} dbl3_t; /* ---------------------------------------------------------------------- */ @@ -87,23 +91,23 @@ PairList::~PairList() void PairList::compute(int eflag, int vflag) { - ev_init(eflag,vflag); + ev_init(eflag, vflag); // get maximum allowed tag. bigint maxtag_one, maxtag; maxtag_one = maxtag = 0; const int nlocal = atom->nlocal; - const tagint * _noalias const tag = atom->tag; + const tagint *_noalias const tag = atom->tag; for (int i = 0; i < nlocal; ++i) maxtag_one = MAX(maxtag_one, tag[i]); MPI_Allreduce(&maxtag_one, &maxtag, 1, MPI_LMP_TAGINT, MPI_MAX, world); const int newton_pair = force->newton_pair; - const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; - dbl3_t * _noalias const f = (dbl3_t *) atom->f[0]; // NOLINT + const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0]; + dbl3_t *_noalias const f = (dbl3_t *) atom->f[0]; // NOLINT - double fpair,epair; - int i,j; + double fpair, epair; + int i, j; int pc = 0; for (int n = 0; n < npairs; ++n) { @@ -113,12 +117,14 @@ void PairList::compute(int eflag, int vflag) if ((par.id1 < 1) || (par.id1 > maxtag)) { if (check_flag) error->all(FLERR, "Invalid pair list atom ID {}", par.id1); - else continue; + else + continue; } if ((par.id2 < 1) || (par.id2 > maxtag)) { if (check_flag) error->all(FLERR, "Invalid pair list atom ID {}", par.id2); - else continue; + else + continue; } i = atom->map(par.id1); @@ -134,14 +140,14 @@ void PairList::compute(int eflag, int vflag) // if id1 is a ghost, we skip if the sum of both ids is even. // if id2 is a ghost, we skip if the sum of both ids is odd. if (newton_pair) { - if ((i >= nlocal) && ((par.id1+par.id2) & 1) == 0) continue; - if ((j >= nlocal) && ((par.id1+par.id2) & 1) == 1) continue; + if ((i >= nlocal) && ((par.id1 + par.id2) & 1) == 0) continue; + if ((j >= nlocal) && ((par.id1 + par.id2) & 1) == 1) continue; } const double dx = x[i].x - x[j].x; const double dy = x[i].y - x[j].y; const double dz = x[i].z - x[j].z; - const double rsq = dx*dx + dy*dy + dz*dz; + const double rsq = dx * dx + dy * dy + dz * dz; fpair = epair = 0.0; if (check_flag) { @@ -150,30 +156,28 @@ void PairList::compute(int eflag, int vflag) } if (rsq < par.cutsq) { - const double r2inv = 1.0/rsq; + const double r2inv = 1.0 / rsq; if (par.style == HARM) { const double r = sqrt(rsq); const double dr = par.param.harm.r0 - r; - fpair = 2.0*par.param.harm.k*dr/r; + fpair = 2.0 * par.param.harm.k * dr / r; - if (eflag_either) - epair = par.param.harm.k*dr*dr - par.offset; + if (eflag_either) epair = par.param.harm.k * dr * dr - par.offset; } else if (par.style == MORSE) { const double r = sqrt(rsq); 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; + fpair = 2.0 * par.param.morse.d0 * par.param.morse.alpha * (dexp * dexp - dexp) / r; - if (eflag_either) - epair = par.param.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; + if (eflag_either) epair = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp) - par.offset; } else if (par.style == LJ126) { - const double r6inv = r2inv*r2inv*r2inv; - const double sig6 = mypow(par.param.lj126.sigma,6); + const double r6inv = r2inv * r2inv * r2inv; + 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) @@ -188,32 +192,31 @@ void PairList::compute(int eflag, int vflag) double r2 = dr * dr; fpair = -par.param.quartic.k / r * (r2 * (ra + rb) + 2.0 * dr * ra * rb); - if (eflag_either) - epair = par.param.quartic.k * r2 * ra * rb; + if (eflag_either) epair = par.param.quartic.k * r2 * ra * rb; } if (newton_pair || i < nlocal) { - f[i].x += dx*fpair; - f[i].y += dy*fpair; - f[i].z += dz*fpair; + f[i].x += dx * fpair; + f[i].y += dy * fpair; + f[i].z += dz * fpair; } if (newton_pair || j < nlocal) { - f[j].x -= dx*fpair; - f[j].y -= dy*fpair; - f[j].z -= dz*fpair; + f[j].x -= dx * fpair; + f[j].y -= dy * fpair; + f[j].z -= dz * fpair; } - if (evflag) ev_tally(i,j,nlocal,newton_pair,epair,0.0,fpair,dx,dy,dz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, epair, 0.0, fpair, dx, dy, dz); } } if (vflag_fdotr) virial_fdotr_compute(); if (check_flag) { int tmp; - MPI_Allreduce(&pc,&tmp,1,MPI_INT,MPI_SUM,world); - if (tmp != 2*npairs) - error->all(FLERR,"Not all pairs processed in pair_style list"); + MPI_Allreduce(&pc, &tmp, 1, MPI_INT, MPI_SUM, world); + if (tmp != 2 * npairs) + error->all(FLERR, "Not all pairs processed in pair_style list: {} vs {}", tmp, 2 * npairs); } } @@ -226,12 +229,11 @@ void PairList::allocate() allocated = 1; int np1 = atom->ntypes + 1; - memory->create(setflag,np1,np1,"pair:setflag"); + memory->create(setflag, np1, np1, "pair:setflag"); for (int i = 1; i < np1; i++) - for (int j = i; j < np1; j++) - setflag[i][j] = 0; + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,np1,np1,"pair:cutsq"); + memory->create(cutsq, np1, np1, "pair:cutsq"); } /* ---------------------------------------------------------------------- @@ -242,16 +244,17 @@ void PairList::settings(int narg, char **arg) { if (narg < 2) utils::missing_cmd_args(FLERR, "pair_style list", error); - cut_global = utils::numeric(FLERR,arg[1],false,lmp); + cut_global = utils::numeric(FLERR, arg[1], false, lmp); int iarg = 2; while (iarg < narg) { - if (strcmp(arg[iarg],"nocheck") == 0) { + if (strcmp(arg[iarg], "nocheck") == 0) { check_flag = 0; ++iarg; - } else if (strcmp(arg[2],"check") == 0) { + } else if (strcmp(arg[2], "check") == 0) { check_flag = 1; ++iarg; - } else error->all(FLERR, "Unknown pair style list keyword: {}", arg[iarg]); + } else + error->all(FLERR, "Unknown pair style list keyword: {}", arg[iarg]); } std::vector mystyles; @@ -260,7 +263,7 @@ void PairList::settings(int narg, char **arg) // read and parse potential file only on MPI rank 0. if (comm->me == 0) { int nharm, nmorse, nlj126, nquartic, nskipped; - FILE *fp = utils::open_potential(arg[0],lmp,nullptr); + FILE *fp = utils::open_potential(arg[0], lmp, nullptr); if (!fp) error->one(FLERR, "Error opening pair list coeffs file {}: {}", arg[0], utils::getsyserror()); TextFileReader reader(fp, "pair list coeffs"); @@ -279,60 +282,62 @@ void PairList::settings(int narg, char **arg) switch (oneparam.style) { - case HARM: - oneparam.param.harm.k = values.next_double(); - oneparam.param.harm.r0 = values.next_double(); - ++nharm; - break; + 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 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 LJ126: + oneparam.param.lj126.epsilon = values.next_double(); + oneparam.param.lj126.sigma = values.next_double(); + ++nlj126; + break; - case QUARTIC: - oneparam.param.quartic.k = values.next_double(); - oneparam.param.quartic.b1 = values.next_double(); - oneparam.param.quartic.b2 = values.next_double(); - if (!values.has_next()) - throw FileReaderException("Must specify individual cutoff for quartic interaction"); - ++nquartic; - break; + case QUARTIC: + oneparam.param.quartic.k = values.next_double(); + oneparam.param.quartic.b1 = values.next_double(); + oneparam.param.quartic.b2 = values.next_double(); + if (!values.has_next()) + throw FileReaderException("Must specify individual cutoff for quartic interaction"); + ++nquartic; + break; - case NONE: // fallthrough - error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}", - utils::trim(line)); - ++nskipped; - break; + case NONE: // fallthrough + error->warning(FLERR, "Skipping unrecognized pair list potential entry: {}", + utils::trim(line)); + ++nskipped; + break; } if (values.has_next()) oneparam.cutsq = square(values.next_double()); else - oneparam.cutsq = cut_global*cut_global; + oneparam.cutsq = cut_global * cut_global; myparams.push_back(oneparam); } } catch (std::exception &e) { - error->one(FLERR,"Error reading pair list coeffs file: {}\n{}", e.what(),line); + error->one(FLERR, "Error reading pair list coeffs file: {}\n{}", e.what(), line); } - utils::logmesg(lmp, "Read {} ({}/{}/{}/{}) interacting pair lines from {}. " - "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, nquartic, arg[0], nskipped); + utils::logmesg(lmp, + "Read {} ({}/{}/{}/{}) interacting pair lines from {}. " + "{} skipped entries.\n", + npairs, nharm, nmorse, nlj126, nquartic, arg[0], nskipped); - memory->create(params,npairs,"pair_list:params"); - memcpy(params, myparams.data(),npairs*sizeof(list_param)); + 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_BYTE, 0, world); + if (comm->me != 0) memory->create(params, npairs, "pair_list:params"); + MPI_Bcast(params, npairs * sizeof(list_param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- @@ -341,22 +346,22 @@ void PairList::settings(int narg, char **arg) void PairList::coeff(int narg, char **arg) { - if (narg < 2) utils::missing_cmd_args(FLERR,"pair_coeff list", error); + if (narg < 2) utils::missing_cmd_args(FLERR, "pair_coeff list", error); 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); + 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); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { setflag[i][j] = 1; count++; } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -365,29 +370,27 @@ void PairList::coeff(int narg, char **arg) void PairList::init_style() { - if (atom->tag_enable == 0) - error->all(FLERR,"Pair style list requires atom IDs"); + if (atom->tag_enable == 0) error->all(FLERR, "Pair style list requires atom IDs"); - if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Pair style list requires an atom map"); + if (atom->map_style == Atom::MAP_NONE) error->all(FLERR, "Pair style list requires an atom map"); if (offset_flag) { - for (int n=0; n < npairs; ++n) { + for (int n = 0; n < npairs; ++n) { list_param &par = params[n]; if (par.style == HARM) { const double dr = sqrt(par.cutsq) - par.param.harm.r0; - par.offset = par.param.harm.k*dr*dr; + par.offset = par.param.harm.k * dr * dr; } 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); + par.offset = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp); } else if (par.style == LJ126) { - const double r6inv = par.cutsq*par.cutsq*par.cutsq; - const double sig6 = mypow(par.param.lj126.sigma,6); - par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); + const double r6inv = par.cutsq * par.cutsq * par.cutsq; + const double sig6 = mypow(par.param.lj126.sigma, 6); + par.offset = 4.0 * par.param.lj126.epsilon * r6inv * (sig6 * sig6 * r6inv - sig6); } else if (par.style == QUARTIC) { // the offset is always 0 at rc @@ -413,10 +416,10 @@ double PairList::init_one(int, int) double PairList::memory_usage() { - double bytes = (double)npairs * sizeof(int); - 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 *)); + double bytes = (double) npairs * sizeof(int); + 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 *)); return bytes; } From 4a06559da5d8ff43f23dbb77fac19b86f0859b18 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 05:18:42 -0400 Subject: [PATCH 16/21] update docs --- doc/src/pair_list.rst | 50 +++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/doc/src/pair_list.rst b/doc/src/pair_list.rst index 76b31623ad..be7ec0c239 100644 --- a/doc/src/pair_list.rst +++ b/doc/src/pair_list.rst @@ -30,15 +30,17 @@ Description """"""""""" Style *list* computes interactions between explicitly listed pairs of -atoms with the option to select functional form and parameters for -each individual pair. Because the parameters are set in the list -file, the pair_coeff command has no parameters (but still needs to be -provided). The *check* and *nocheck* keywords enable/disable a test -that checks whether all listed bonds were present and computed. +atoms with the option to select functional form and parameters for each +individual pair. Because the parameters are set in the list file, the +pair_coeff command has no parameters (but still needs to be provided). +The *check* and *nocheck* keywords enable/disable tests that checks +whether all listed pairs of atom IDs were present and the interactions +computed. If *nocheck* is set and either atom ID is not present, the +interaction is skipped. This pair style can be thought of as a hybrid between bonded, -non-bonded, and restraint interactions. It will typically be used as -an additional interaction within the *hybrid/overlay* pair style. It +non-bonded, and restraint interactions. It will typically be used as an +additional interaction within the *hybrid/overlay* pair style. It currently supports three interaction styles: a 12-6 Lennard-Jones, a Morse and a harmonic potential. @@ -55,10 +57,10 @@ The format of the list file is as follows: ID2 = atom ID of second atom style = style of interaction coeffs = list of coeffs - cutoff = cutoff for interaction (optional) + cutoff = cutoff for interaction (optional, except for style *quartic*) -The cutoff parameter is optional. If not specified, the global cutoff -is used. +The cutoff parameter is optional for all but the *quartic* interactions. +If it is not specified, the global cutoff is used. Here is an example file: @@ -120,7 +122,7 @@ and the coefficients: * :math:`b_2` (distance units) * :math:`r_c` (distance units) -Note that the global cutoff specified by *pair_style list* command is ignored, unless :math:`r_c` is not specified. In this case, :math:`r_c` equals the sqare root of the globle cutoff. +Note that the per list entry cutoff :math:`r_c` is **required** for *quartic* interactions. ---------- @@ -136,8 +138,9 @@ pair style. The :doc:`pair_modify ` table and tail options are not relevant for this pair style. -This pair style does not write its information to :doc:`binary restart files `, so pair_style and pair_coeff commands need -to be specified in an input script that reads a restart file. +This pair style does not write its information to :doc:`binary restart +files `, so pair_style and pair_coeff commands need to be +specified in an input script that reads a restart file. This pair style can only be used via the *pair* keyword of the :doc:`run_style respa ` command. It does not support the @@ -149,17 +152,18 @@ Restrictions """""""""""" This pair style does not use a neighbor list and instead identifies -atoms by their IDs. This has two consequences: 1) The cutoff has to be -chosen sufficiently large, so that the second atom of a pair has to be -a ghost atom on the same node on which the first atom is local; -otherwise the interaction will be skipped. You can use the *check* -option to detect, if interactions are missing. 2) Unlike other pair -styles in LAMMPS, an atom I will not interact with multiple images of -atom J (assuming the images are within the cutoff distance), but only -with the nearest image. +atoms by their IDs. This has two consequences: 1) The cutoff has to be +chosen sufficiently large, so that the second atom of a pair has to be a +ghost atom on the same node on which the first atom is local; otherwise +the interaction will be skipped. You can use the *check* option to +detect, if interactions are missing. 2) Unlike other pair styles in +LAMMPS, an atom I will not interact with multiple images of atom J +(assuming the images are within the cutoff distance), but only with the +closest image. -This style is part of the MISC package. It is only enabled if -LAMMPS is build with that package. See the :doc:`Build package ` page on for more info. +This style is part of the MISC package. It is only enabled if LAMMPS is +build with that package. See the :doc:`Build package ` +page on for more info. Related commands """""""""""""""" From f2b3d8e8c35e18ef605922a8956bca5e0326676f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 07:02:44 -0400 Subject: [PATCH 17/21] use morse interaction energy from bond style morse not pair style morse --- doc/src/pair_list.rst | 5 +++-- src/MISC/pair_list.cpp | 8 ++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/doc/src/pair_list.rst b/doc/src/pair_list.rst index be7ec0c239..762afb888c 100644 --- a/doc/src/pair_list.rst +++ b/doc/src/pair_list.rst @@ -88,7 +88,7 @@ The style *morse* computes pairwise interactions with the formula .. math:: - E = D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right] \qquad r < r_c + E = D_0 \left[ 1 - e^{-\alpha (r - r_0)} \right]^2 \qquad r < r_c and the coefficients: @@ -171,8 +171,9 @@ Related commands :doc:`pair_coeff `, :doc:`pair_style hybrid/overlay `, :doc:`pair_style lj/cut `, -:doc:`pair_style morse `, +:doc:`bond_style morse `, :doc:`bond_style harmonic ` +:doc:`bond_style quartic ` Default """"""" diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index e826ff0693..f6a27538ab 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -168,11 +168,11 @@ void PairList::compute(int eflag, int vflag) } else if (par.style == MORSE) { const double r = sqrt(rsq); - const double dr = par.param.morse.r0 - r; - const double dexp = exp(par.param.morse.alpha * dr); + const double dr = r - par.param.morse.r0; + 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.param.morse.d0 * (dexp * dexp - 2.0 * dexp) - par.offset; + if (eflag_either) epair = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp + 1.0) - par.offset; } else if (par.style == LJ126) { @@ -385,7 +385,7 @@ void PairList::init_style() } 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); + par.offset = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp - 1.0); } else if (par.style == LJ126) { const double r6inv = par.cutsq * par.cutsq * par.cutsq; From 3ff203b70521ee995eb529db17549035aa689450 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 07:03:13 -0400 Subject: [PATCH 18/21] add unit test for pair style list comparing it to pair lj/cut + harmonic/morse bonds --- unittest/force-styles/CMakeLists.txt | 4 + unittest/force-styles/test_pair_list.cpp | 99 ++++++++++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 unittest/force-styles/test_pair_list.cpp diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt index 294ece3a1f..878492d6c7 100644 --- a/unittest/force-styles/CMakeLists.txt +++ b/unittest/force-styles/CMakeLists.txt @@ -260,3 +260,7 @@ if(MLIAP_ENABLE_PYTHON AND (NOT WIN32)) set_tests_properties(TestMliapPyUnified PROPERTIES ENVIRONMENT "PYTHONPATH=${LAMMPS_PYTHON_DIR};PYTHONDONTWRITEBYTECODE=1") endif() +add_executable(test_pair_list test_pair_list.cpp) +target_link_libraries(test_pair_list PRIVATE lammps GTest::GMockMain) +add_test(NAME TestPairList COMMAND test_pair_list) + diff --git a/unittest/force-styles/test_pair_list.cpp b/unittest/force-styles/test_pair_list.cpp new file mode 100644 index 0000000000..32a38d8708 --- /dev/null +++ b/unittest/force-styles/test_pair_list.cpp @@ -0,0 +1,99 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "library.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +const char parms[] = "print \"\"\"\n" + "1 2 lj126 0.3 3.5 4.0\n" + "2 3 lj126 0.3 3.5 5.0\n" + "1 4 harmonic 10.0 7.0\n" + "2 4 harmonic 10.0 7.0\n" + "3 4 morse 10.0 0.3 6.5\n" + "\"\"\" file list.param\n"; + +const char first[] = "units real\n" + "atom_style bond\n" + "atom_modify map array\n" + "boundary f f f\n" + "special_bonds lj/coul 0.0 1.0 1.0\n" + "region box block -5 5 -5 5 -5 5\n" + "create_box 1 box bond/types 2 extra/bond/per/atom 4\n" + "create_atoms 1 single -2.0 0.0 0.0\n" + "create_atoms 1 single 0.0 1.0 0.0\n" + "create_atoms 1 single 4.0 1.0 0.0\n" + "create_atoms 1 single 4.0 -4.0 -4.0\n" + "create_bonds single/bond 1 1 4\n" + "create_bonds single/bond 1 2 4\n" + "create_bonds single/bond 2 3 4\n" + "mass 1 10.0\n" + "velocity all create 10.0 87287 loop geom\n"; + +const char second[] = "timestep 0.2\n" + "fix 1 all nve\n" + "run 2 post no\n"; + +static constexpr double EPSILON = 1.0e-10; + +namespace LAMMPS_NS { + +TEST(PairList, ListVsPairBond) +{ + if (!lammps_config_has_package("MOLECULE")) GTEST_SKIP(); + + const char *lmpargv[] = {"melt", "-log", "none", "-nocite"}; + int lmpargc = sizeof(lmpargv) / sizeof(const char *); + + ::testing::internal::CaptureStdout(); + void *ljmelt = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr); + lmpargv[0] = "plist"; + void *plist = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr); + + lammps_commands_string(ljmelt, first); + lammps_command(ljmelt, "pair_style lj/cut 5.0"); + lammps_command(ljmelt, "pair_coeff * * 0.3 3.5"); + lammps_command(ljmelt, "bond_style hybrid harmonic morse"); + lammps_command(ljmelt, "bond_coeff 1 harmonic 10.0 7.0"); + lammps_command(ljmelt, "bond_coeff 2 morse 10.0 0.3 6.5"); + + lammps_commands_string(ljmelt, second); + + lammps_command(plist, parms); + lammps_commands_string(plist, first); + lammps_command(plist, "pair_style list list.param 10.0"); + lammps_command(plist, "pair_coeff * *"); + lammps_command(plist, "bond_style zero"); + lammps_command(plist, "bond_coeff * 2.0"); + lammps_commands_string(plist, second); + ::testing::internal::GetCapturedStdout(); + + double lj_pe = lammps_get_thermo(ljmelt, "pe"); + double ml_pe = lammps_get_thermo(plist, "pe"); + EXPECT_NEAR(lj_pe, ml_pe, EPSILON); + double lj_ke = lammps_get_thermo(ljmelt, "ke"); + double ml_ke = lammps_get_thermo(plist, "ke"); + EXPECT_NEAR(lj_ke, ml_ke, EPSILON); + double lj_press = lammps_get_thermo(ljmelt, "press"); + double ml_press = lammps_get_thermo(plist, "press"); + EXPECT_NEAR(lj_press, ml_press, EPSILON); + + ::testing::internal::CaptureStdout(); + lammps_command(plist, "shell rm list.param"); + lammps_close(ljmelt); + lammps_close(plist); + ::testing::internal::GetCapturedStdout(); +} + +} // namespace LAMMPS_NS From 7f90e43d0a31f8ad403d3d5712a760b0e93a158d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 14:34:52 -0400 Subject: [PATCH 19/21] must skip test for pair style list, if MISC package is not available --- unittest/force-styles/test_pair_list.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/unittest/force-styles/test_pair_list.cpp b/unittest/force-styles/test_pair_list.cpp index 32a38d8708..9e249613f3 100644 --- a/unittest/force-styles/test_pair_list.cpp +++ b/unittest/force-styles/test_pair_list.cpp @@ -52,6 +52,7 @@ namespace LAMMPS_NS { TEST(PairList, ListVsPairBond) { if (!lammps_config_has_package("MOLECULE")) GTEST_SKIP(); + if (!lammps_config_has_package("MISC")) GTEST_SKIP(); const char *lmpargv[] = {"melt", "-log", "none", "-nocite"}; int lmpargc = sizeof(lmpargv) / sizeof(const char *); From 82d6c1de9927ebb481f15256b8f46fe8317445c3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Oct 2022 15:48:01 -0400 Subject: [PATCH 20/21] recover MDI library build with windows cross-compilation --- cmake/Modules/Packages/MDI.cmake | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/cmake/Modules/Packages/MDI.cmake b/cmake/Modules/Packages/MDI.cmake index 82b036a583..c1e25347af 100644 --- a/cmake/Modules/Packages/MDI.cmake +++ b/cmake/Modules/Packages/MDI.cmake @@ -49,6 +49,14 @@ if(DOWNLOAD_MDI) set(MDI_USE_PYTHON_PLUGINS ON) endif() endif() + # python plugins are not supported and thus must be always off on Windows + if(CMAKE_SYSTEM_NAME STREQUAL "Windows") + unset(Python_Development_FOUND) + set(MDI_USE_PYTHON_PLUGINS OFF) + if(CMAKE_CROSSCOMPILING) + set(CMAKE_INSTALL_LIBDIR lib) + endif() + endif() # download/ build MDI library # always build static library with -fpic @@ -57,8 +65,9 @@ if(DOWNLOAD_MDI) ExternalProject_Add(mdi_build URL ${MDI_URL} URL_MD5 ${MDI_MD5} + PREFIX ${CMAKE_CURRENT_BINARY_DIR}/mdi_build_ext CMAKE_ARGS - -DCMAKE_INSTALL_PREFIX= + -DCMAKE_INSTALL_PREFIX=${CMAKE_CURRENT_BINARY_DIR}/mdi_build_ext -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM} @@ -70,21 +79,22 @@ if(DOWNLOAD_MDI) -Dplugins=ON -Dpython_plugins=${MDI_USE_PYTHON_PLUGINS} UPDATE_COMMAND "" - BUILD_BYPRODUCTS "/${CMAKE_INSTALL_LIBDIR}/mdi/${CMAKE_STATIC_LIBRARY_PREFIX}mdi${CMAKE_STATIC_LIBRARY_SUFFIX}" + INSTALL_COMMAND ${CMAKE_COMMAND} --build ${CMAKE_CURRENT_BINARY_DIR}/mdi_build_ext/src/mdi_build-build --target install + BUILD_BYPRODUCTS "${CMAKE_CURRENT_BINARY_DIR}/mdi_build_ext/${CMAKE_INSTALL_LIBDIR}/mdi/${CMAKE_STATIC_LIBRARY_PREFIX}mdi${CMAKE_STATIC_LIBRARY_SUFFIX}" ) # where is the compiled library? - ExternalProject_get_property(mdi_build INSTALL_DIR) + ExternalProject_get_property(mdi_build PREFIX) # workaround for older CMake versions - file(MAKE_DIRECTORY ${INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/mdi) - file(MAKE_DIRECTORY ${INSTALL_DIR}/include/mdi) + file(MAKE_DIRECTORY ${PREFIX}/${CMAKE_INSTALL_LIBDIR}/mdi) + file(MAKE_DIRECTORY ${PREFIX}/include/mdi) # create imported target for the MDI library add_library(LAMMPS::MDI UNKNOWN IMPORTED) add_dependencies(LAMMPS::MDI mdi_build) set_target_properties(LAMMPS::MDI PROPERTIES - IMPORTED_LOCATION "${INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/mdi/${CMAKE_STATIC_LIBRARY_PREFIX}mdi${CMAKE_STATIC_LIBRARY_SUFFIX}" - INTERFACE_INCLUDE_DIRECTORIES ${INSTALL_DIR}/include/mdi + IMPORTED_LOCATION "${PREFIX}/${CMAKE_INSTALL_LIBDIR}/mdi/${CMAKE_STATIC_LIBRARY_PREFIX}mdi${CMAKE_STATIC_LIBRARY_SUFFIX}" + INTERFACE_INCLUDE_DIRECTORIES ${PREFIX}/include/mdi ) set(MDI_DEP_LIBS "") From 80da4c307c5a18026bfd47549255825cbb684781 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 5 Oct 2022 10:52:17 -0400 Subject: [PATCH 21/21] silence compiler warnings, avoid integer or buffer overflows --- fortran/lammps.f90 | 12 +++---- .../fortran/test_fortran_extract_variable.f90 | 36 +++++++++---------- .../fortran/test_fortran_gather_scatter.f90 | 1 - unittest/fortran/wrap_extract_variable.cpp | 4 +-- 4 files changed, 23 insertions(+), 30 deletions(-) diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index d37b12c7db..880eef9b67 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -1134,7 +1134,7 @@ CONTAINS TYPE(lammps_variable_data) :: variable_data TYPE(c_ptr) :: Cptr, Cname, Cgroup, Cveclength - INTEGER :: length, i + INTEGER(c_size_t) :: length, i CHARACTER(KIND=c_char, LEN=1), DIMENSION(:), POINTER :: Cstring INTEGER(c_int) :: datatype REAL(c_double), POINTER :: double => NULL() @@ -1370,7 +1370,6 @@ CONTAINS INTEGER(c_int) :: ndata TYPE(c_ptr) :: Cdata, Cname, Cids INTEGER(c_int), PARAMETER :: Ctype = 0_c_int - CHARACTER(LEN=100) :: error_msg IF (count /= 1 .AND. count /= 3) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & @@ -1401,7 +1400,6 @@ CONTAINS INTEGER(c_int) :: ndata TYPE(c_ptr) :: Cdata, Cname, Cids INTEGER(c_int), PARAMETER :: Ctype = 1_c_int - CHARACTER(LEN=100) :: error_msg IF (count /= 1 .AND. count /= 3) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & @@ -1494,7 +1492,6 @@ CONTAINS INTEGER(c_int), PARAMETER :: Ctype = 0_c_int INTEGER(c_int) :: Cndata, Ccount TYPE(c_ptr) :: Cdata, Cname, Cids - CHARACTER(LEN=100) :: error_msg Cndata = SIZE(ids, KIND=c_int) Ccount = SIZE(data, KIND=c_int) / Cndata @@ -1519,7 +1516,6 @@ CONTAINS INTEGER(c_int), PARAMETER :: Ctype = 1_c_int INTEGER(c_int) :: Cndata, Ccount TYPE(c_ptr) :: Cdata, Cname, Cids - CHARACTER(LEN=100) :: error_msg Cndata = SIZE(ids, KIND=c_int) Ccount = SIZE(data, KIND=c_int) / Cndata @@ -1628,7 +1624,7 @@ CONTAINS INTEGER(c_int) :: Cidx, Csuccess TYPE(c_ptr) :: Cptr CHARACTER(LEN=1,KIND=c_char), TARGET :: Cbuffer(LEN(buffer)+1) - INTEGER :: i, strlen + INTEGER(c_size_t) :: i, strlen Cidx = idx - 1 Cptr = C_LOC(Cbuffer(1)) @@ -1698,8 +1694,8 @@ CONTAINS CLASS(lammps), INTENT(IN) :: self CHARACTER(LEN=*), INTENT(OUT) :: buffer INTEGER, INTENT(OUT), OPTIONAL :: status - INTEGER(c_int) :: buflen, Cstatus, i - INTEGER(c_size_t) :: length + INTEGER(c_int) :: buflen, Cstatus + INTEGER(c_size_t) :: i, length TYPE(c_ptr) :: Cptr CHARACTER(LEN=1, KIND=c_char), POINTER :: c_string(:) diff --git a/unittest/fortran/test_fortran_extract_variable.f90 b/unittest/fortran/test_fortran_extract_variable.f90 index ded3743409..e373228a0f 100644 --- a/unittest/fortran/test_fortran_extract_variable.f90 +++ b/unittest/fortran/test_fortran_extract_variable.f90 @@ -33,8 +33,7 @@ CONTAINS CHARACTER(LEN=256) :: test_input_directory TYPE(c_ptr) :: c_test_input_directory, c_absolute_path, c_filename CHARACTER(LEN=1,KIND=c_char), DIMENSION(:), POINTER :: F_absolute_path - INTEGER :: i - INTEGER(c_size_t) :: length + INTEGER(c_size_t) :: i, length test_input_directory = lmp%extract_variable('input_dir') c_test_input_directory = f2c_string(test_input_directory) @@ -91,10 +90,10 @@ FUNCTION f_lammps_with_C_args(argc, argv) BIND(C) TYPE(c_ptr) :: f_lammps_with_C_args CHARACTER(LEN=ARG_LENGTH), DIMENSION(argc) :: args CHARACTER(LEN=1,KIND=c_char), DIMENSION(:), POINTER :: Cstr - INTEGER :: i, length, j + INTEGER(c_size_t):: i, length, j INTERFACE - FUNCTION c_strlen (str) BIND(C,name='strlen') + FUNCTION c_strlen(str) BIND(C,name='strlen') IMPORT :: c_ptr, c_size_t IMPLICIT NONE TYPE(c_ptr), INTENT(IN), VALUE :: str @@ -126,7 +125,7 @@ SUBROUTINE f_lammps_close() BIND(C) lmp%handle = c_null_ptr END SUBROUTINE f_lammps_close -SUBROUTINE f_lammps_setup_extract_variable () BIND(C) +SUBROUTINE f_lammps_setup_extract_variable() BIND(C) USE LIBLAMMPS USE keepstuff, ONLY : lmp, big_input, cont_input, more_input, pair_input USE keepvar, ONLY : absolute_path @@ -173,7 +172,7 @@ SUBROUTINE f_lammps_setup_extract_variable () BIND(C) CALL lmp%command("run 0") ! so c_COM and v_center have values END SUBROUTINE f_lammps_setup_extract_variable -FUNCTION f_lammps_extract_variable_index_1 () BIND(C) +FUNCTION f_lammps_extract_variable_index_1() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -189,7 +188,7 @@ FUNCTION f_lammps_extract_variable_index_1 () BIND(C) END IF END FUNCTION f_lammps_extract_variable_index_1 -FUNCTION f_lammps_extract_variable_index_2 () BIND(C) +FUNCTION f_lammps_extract_variable_index_2() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -205,7 +204,7 @@ FUNCTION f_lammps_extract_variable_index_2 () BIND(C) END IF END FUNCTION f_lammps_extract_variable_index_2 -FUNCTION f_lammps_extract_variable_loop () BIND(C) +FUNCTION f_lammps_extract_variable_loop() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -217,7 +216,7 @@ FUNCTION f_lammps_extract_variable_loop () BIND(C) READ(loop,*) f_lammps_extract_variable_loop END FUNCTION f_lammps_extract_variable_loop -FUNCTION f_lammps_extract_variable_loop_pad () BIND(C) +FUNCTION f_lammps_extract_variable_loop_pad() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -230,7 +229,7 @@ FUNCTION f_lammps_extract_variable_loop_pad () BIND(C) f_lammps_extract_variable_loop_pad = f2c_string(loop) END FUNCTION f_lammps_extract_variable_loop_pad -FUNCTION f_lammps_extract_variable_world () BIND(C) +FUNCTION f_lammps_extract_variable_world() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -243,7 +242,7 @@ FUNCTION f_lammps_extract_variable_world () BIND(C) f_lammps_extract_variable_world = f2c_string(world) END FUNCTION f_lammps_extract_variable_world -FUNCTION f_lammps_extract_variable_universe () BIND(C) +FUNCTION f_lammps_extract_variable_universe() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -256,7 +255,7 @@ FUNCTION f_lammps_extract_variable_universe () BIND(C) f_lammps_extract_variable_universe = f2c_string(universe) END FUNCTION f_lammps_extract_variable_universe -FUNCTION f_lammps_extract_variable_uloop () BIND(C) +FUNCTION f_lammps_extract_variable_uloop() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -268,7 +267,7 @@ FUNCTION f_lammps_extract_variable_uloop () BIND(C) READ(uloop,*) f_lammps_extract_variable_uloop END FUNCTION f_lammps_extract_variable_uloop -FUNCTION f_lammps_extract_variable_string () BIND(C) +FUNCTION f_lammps_extract_variable_string() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -281,7 +280,7 @@ FUNCTION f_lammps_extract_variable_string () BIND(C) f_lammps_extract_variable_string = f2c_string(string) END FUNCTION f_lammps_extract_variable_string -FUNCTION f_lammps_extract_variable_format () BIND(C) +FUNCTION f_lammps_extract_variable_format() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -294,7 +293,7 @@ FUNCTION f_lammps_extract_variable_format () BIND(C) f_lammps_extract_variable_format = f2c_string(form) END FUNCTION f_lammps_extract_variable_format -FUNCTION f_lammps_extract_variable_format_pad () BIND(C) +FUNCTION f_lammps_extract_variable_format_pad() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -307,7 +306,7 @@ FUNCTION f_lammps_extract_variable_format_pad () BIND(C) f_lammps_extract_variable_format_pad = f2c_string(form) END FUNCTION f_lammps_extract_variable_format_pad -FUNCTION f_lammps_extract_variable_getenv () BIND(C) +FUNCTION f_lammps_extract_variable_getenv() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -320,7 +319,7 @@ FUNCTION f_lammps_extract_variable_getenv () BIND(C) f_lammps_extract_variable_getenv = f2c_string(string) END FUNCTION f_lammps_extract_variable_getenv -FUNCTION f_lammps_extract_variable_file () BIND(C) +FUNCTION f_lammps_extract_variable_file() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr USE LIBLAMMPS USE keepstuff, ONLY : lmp @@ -346,12 +345,11 @@ FUNCTION f_lammps_extract_variable_atomfile(i) BIND(C) f_lammps_extract_variable_atomfile = atom_data(i) END FUNCTION f_lammps_extract_variable_atomfile -FUNCTION f_lammps_extract_variable_python(i) BIND(C) +FUNCTION f_lammps_extract_variable_python() BIND(C) USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double USE LIBLAMMPS USE keepstuff, ONLY : lmp IMPLICIT NONE - INTEGER(c_int), INTENT(IN), VALUE :: i REAL(c_double) :: f_lammps_extract_variable_python f_lammps_extract_variable_python = lmp%extract_variable('py') diff --git a/unittest/fortran/test_fortran_gather_scatter.f90 b/unittest/fortran/test_fortran_gather_scatter.f90 index 69bb0e030f..ec1880c908 100644 --- a/unittest/fortran/test_fortran_gather_scatter.f90 +++ b/unittest/fortran/test_fortran_gather_scatter.f90 @@ -190,7 +190,6 @@ SUBROUTINE f_lammps_scatter_atoms_subset_mask() BIND(C) INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: all_masks INTEGER(c_int), DIMENSION(*), PARAMETER :: tags = [3,1] INTEGER(c_int), DIMENSION(2) :: masks - INTEGER(c_int) :: swap CALL lmp%gather_atoms('mask', 1_c_int, all_masks) diff --git a/unittest/fortran/wrap_extract_variable.cpp b/unittest/fortran/wrap_extract_variable.cpp index 0c1ffcc37e..1082a381bb 100644 --- a/unittest/fortran/wrap_extract_variable.cpp +++ b/unittest/fortran/wrap_extract_variable.cpp @@ -165,7 +165,7 @@ TEST_F(LAMMPS_extract_variable, format) { f_lammps_setup_extract_variable(); int i; - char str[10]; + char str[16]; char *fstr; for (i = 1; i <= 10; i++) { std::sprintf(str, "%.6G", std::exp(i)); @@ -180,7 +180,7 @@ TEST_F(LAMMPS_extract_variable, format_pad) { f_lammps_setup_extract_variable(); int i; - char str[10]; + char str[16]; char *fstr; for (i = 1; i <= 10; i++) { std::sprintf(str, "%08.6G", std::exp(i));