diff --git a/doc/src/pair_list.rst b/doc/src/pair_list.rst index 321bf61be9..762afb888c 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: @@ -69,6 +71,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 @@ -85,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: @@ -106,6 +109,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 per list entry cutoff :math:`r_c` is **required** for *quartic* interactions. + ---------- Mixing, shift, table, tail correction, restart, rRESPA info @@ -120,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 @@ -133,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 """""""""""""""" @@ -151,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 f6ac2e3190..f6a27538ab 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 @@ -22,38 +21,46 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "math_special.h" #include "memory.h" #include "text_file_reader.h" -#include #include +#include #include #include using namespace LAMMPS_NS; +using MathSpecial::square; -enum { NONE = 0, HARM, MORSE, LJ126 }; +enum { NONE = 0, HARM, MORSE, LJ126, QUARTIC }; +// clang-format off static std::map stylename = { - { "none", NONE }, - { "harmonic", HARM }, - { "morse", MORSE }, - { "lj126", LJ126 } + {"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; /* ---------------------------------------------------------------------- */ @@ -84,19 +91,42 @@ 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 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 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); - double fpair,epair; - int i,j; + 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 + + double fpair, epair; + 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); @@ -110,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) { @@ -126,61 +156,67 @@ 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; + 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) { - 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; + 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) - 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; + 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; - 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); } } @@ -191,14 +227,13 @@ 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++) - setflag[i][j] = 0; + 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"); } /* ---------------------------------------------------------------------- @@ -207,13 +242,19 @@ 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; + cut_global = utils::numeric(FLERR, arg[1], false, lmp); + 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; @@ -221,13 +262,15 @@ 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; - FILE *fp = utils::open_potential(arg[0],lmp,nullptr); - TextFileReader reader(fp,"pair list coeffs"); - npairs = nharm = nmorse = nlj126 = nskipped = 0; + int nharm, nmorse, nlj126, nquartic, nskipped; + 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"); + npairs = nharm = nmorse = nlj126 = nquartic = nskipped = 0; + char *line; try { - char *line; while ((line = reader.next_line())) { ValueTokenizer values(line); list_param oneparam; @@ -239,75 +282,86 @@ 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 NONE: // fallthrough - error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}", - utils::trim(line)); - ++nskipped; - 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; } if (values.has_next()) - oneparam.cutsq = values.next_double(); + 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: {}", 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, 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); } /* ---------------------------------------------------------------------- - 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; - 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"); } /* ---------------------------------------------------------------------- @@ -316,29 +370,31 @@ 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 - 1.0); } 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 + par.offset = 0.0; } } } @@ -360,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; } diff --git a/src/MISC/pair_list.h b/src/MISC/pair_list.h index 2828a08547..dab8ba3fc9 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; + }; + union param_u { harm_p harm; morse_p morse; lj126_p lj126; + quartic_p quartic; }; struct list_param { 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