diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 514785c15c..0629afbf08 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -26,6 +26,7 @@ OPT. * :doc:`none ` * :doc:`zero ` * :doc:`hybrid (k) ` + * :doc:`hybrid/molecular ` * :doc:`hybrid/overlay (k) ` * :doc:`hybrid/scaled ` * :doc:`kim ` @@ -34,7 +35,6 @@ OPT. * * * - * * :doc:`adp (ko) ` * :doc:`agni (o) ` * :doc:`aip/water/2dm (t) ` diff --git a/doc/src/pair_hybrid.rst b/doc/src/pair_hybrid.rst index a33991e43e..4108aa6154 100644 --- a/doc/src/pair_hybrid.rst +++ b/doc/src/pair_hybrid.rst @@ -1,5 +1,6 @@ .. index:: pair_style hybrid .. index:: pair_style hybrid/kk +.. index:: pair_style hybrid/molecular .. index:: pair_style hybrid/overlay .. index:: pair_style hybrid/overlay/kk .. index:: pair_style hybrid/scaled @@ -9,6 +10,9 @@ pair_style hybrid command Accelerator Variants: *hybrid/kk* +pair_style hybrid/molecular command +=================================== + pair_style hybrid/overlay command ================================= @@ -23,6 +27,7 @@ Syntax .. code-block:: LAMMPS pair_style hybrid style1 args style2 args ... + pair_style hybrid/molecular factor1 style1 args factor2 style 2 args pair_style hybrid/overlay style1 args style2 args ... pair_style hybrid/scaled factor1 style1 args factor2 style 2 args ... @@ -47,6 +52,10 @@ Examples pair_coeff * * tersoff Si.tersoff Si pair_coeff * * sw Si.sw Si + pair_style hybrid/molecular lj/cut 2.5 lj/cut 2.5 + pair_coeff * * lj/cut 1 1.0 1.0 + pair_coeff * * lj/cut 2 1.5 1.0 + variable one equal ramp(1.0,0.0) variable two equal 1.0-v_one pair_style hybrid/scaled v_one lj/cut 2.5 v_two morse 2.5 @@ -56,17 +65,26 @@ Examples Description """"""""""" -The *hybrid*, *hybrid/overlay*, and *hybrid/scaled* styles enable the -use of multiple pair styles in one simulation. With the *hybrid* style, -exactly one pair style is assigned to each pair of atom types. With the -*hybrid/overlay* and *hybrid/scaled* styles, one or more pair styles can -be assigned to each pair of atom types. The assignment of pair styles -to type pairs is made via the :doc:`pair_coeff ` command. -The major difference between the *hybrid/overlay* and *hybrid/scaled* -styles is that the *hybrid/scaled* adds a scale factor for each -sub-style contribution to forces, energies and stresses. Because of the -added complexity, the *hybrid/scaled* style has more overhead and thus -may be slower than *hybrid/overlay*. +The *hybrid*, *hybrid/overlay*, *hybrid/molecular*, and *hybrid/scaled* +styles enable the use of multiple pair styles in one simulation. With +the *hybrid* style, exactly one pair style is assigned to each pair of +atom types. With the *hybrid/overlay* and *hybrid/scaled* styles, one +or more pair styles can be assigned to each pair of atom types. With +the hybrid/molecular style, pair styles are assigned to either intra- +or inter-molecular interactions. + +The assignment of pair styles to type pairs is made via the +:doc:`pair_coeff ` command. The major difference between +the *hybrid/overlay* and *hybrid/scaled* styles is that the +*hybrid/scaled* adds a scale factor for each sub-style contribution to +forces, energies and stresses. Because of the added complexity, the +*hybrid/scaled* style has more overhead and thus may be slower than +*hybrid/overlay*. + +The *hybrid/molecular* pair style accepts *only* two sub-styles: the +first is assigned to intra-molecular interactions (i.e. both atoms +have the same molecule ID), the second to inter-molecular interactions +(i.e. interacting atoms have different molecule IDs). Here are two examples of hybrid simulations. The *hybrid* style could be used for a simulation of a metal droplet on a LJ surface. The metal @@ -476,6 +494,8 @@ the same or else LAMMPS will generate an error. Pair style *hybrid/scaled* currently only works for non-accelerated pair styles and pair styles from the OPT package. +Pair style *hybrid/molecular* is not compatible with manybody potentials. + When using pair styles from the GPU package they must not be listed multiple times. LAMMPS will detect this and abort. diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 74dfce6b01..13dd922043 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -108,6 +108,7 @@ accelerated styles exist. * :doc:`none ` - turn off pairwise interactions * :doc:`hybrid ` - multiple styles of pairwise interactions +* :doc:`hybrid/molecular ` - different pair styles for intra- and inter-molecular interactions * :doc:`hybrid/overlay ` - multiple styles of superposed pairwise interactions * :doc:`hybrid/scaled ` - multiple styles of scaled superposed pairwise interactions * :doc:`zero ` - neighbor list but no interactions diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp index 4bdd58eead..878ee98917 100644 --- a/src/neigh_list.cpp +++ b/src/neigh_list.cpp @@ -148,6 +148,7 @@ void NeighList::post_constructor(NeighRequest *nq) copy = nq->copy; trim = nq->trim; id = nq->id; + molskip = nq->molskip; if (nq->copy) { listcopy = neighbor->lists[nq->copylist]; @@ -157,14 +158,16 @@ void NeighList::post_constructor(NeighRequest *nq) if (nq->skip) { listskip = neighbor->lists[nq->skiplist]; - int ntypes = atom->ntypes; - iskip = new int[ntypes+1]; - memory->create(ijskip,ntypes+1,ntypes+1,"neigh_list:ijskip"); - int i,j; - for (i = 1; i <= ntypes; i++) iskip[i] = nq->iskip[i]; - for (i = 1; i <= ntypes; i++) - for (j = 1; j <= ntypes; j++) - ijskip[i][j] = nq->ijskip[i][j]; + if (!molskip) { + int ntypes = atom->ntypes; + iskip = new int[ntypes+1]; + memory->create(ijskip,ntypes+1,ntypes+1,"neigh_list:ijskip"); + int i,j; + for (i = 1; i <= ntypes; i++) iskip[i] = nq->iskip[i]; + for (i = 1; i <= ntypes; i++) + for (j = 1; j <= ntypes; j++) + ijskip[i][j] = nq->ijskip[i][j]; + } } if (nq->halffull) diff --git a/src/neigh_list.h b/src/neigh_list.h index 01e9a1c6a7..4c592772ad 100644 --- a/src/neigh_list.h +++ b/src/neigh_list.h @@ -46,6 +46,7 @@ class NeighList : protected Pointers { int kk2cpu; // 1 if this list is copied from Kokkos to CPU int copymode; // 1 if this is a Kokkos on-device copy int id; // copied from neighbor list request + int molskip; // 1/2 if this is an intra-/inter-molecular skip list // data structs to store neighbor pairs I,J and associated values diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp index 280b5c54ae..262bc5254f 100644 --- a/src/neigh_request.cpp +++ b/src/neigh_request.cpp @@ -76,6 +76,7 @@ NeighRequest::NeighRequest(LAMMPS *_lmp) : Pointers(_lmp) skip = 0; iskip = nullptr; ijskip = nullptr; + molskip = 0; // only set when command = 1; @@ -183,6 +184,8 @@ int NeighRequest::identical(NeighRequest *other) int NeighRequest::same_skip(NeighRequest *other) { + if (molskip != other->molskip) return 0; + const int ntypes = atom->ntypes; int same = 1; @@ -307,6 +310,12 @@ void NeighRequest::set_skip(int *_iskip, int **_ijskip) ijskip = _ijskip; } +void NeighRequest::set_molskip(int _molskip) +{ + skip = 1; + molskip = _molskip; +} + void NeighRequest::enable_full() { half = 0; diff --git a/src/neigh_request.h b/src/neigh_request.h index fa57922c93..caa9e05ad6 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -29,6 +29,9 @@ class NeighRequest : protected Pointers { friend class NPairSkipTrimIntel; friend class FixIntel; + public: + enum { REGULAR, INTRA, INTER }; + protected: void *requestor; // class that made request int requestor_instance; // instance of that class (only Fix, Compute, Pair) @@ -88,6 +91,7 @@ class NeighRequest : protected Pointers { int skip; // 1 if this list skips atom types from another list int *iskip; // iskip[i] if atoms of type I are not in list int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list + int molskip; // 0 reqular list, 1 keep only intra-molecular entries, 2 keep inter-molecular // command_style only set if command = 1 // allows print_pair_info() to access command name @@ -137,6 +141,7 @@ class NeighRequest : protected Pointers { void set_kokkos_device(int); void set_kokkos_host(int); void set_skip(int *, int **); + void set_molskip(int); void enable_full(); void enable_ghost(); void enable_intel(); diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 63d14acb9a..59ddf87698 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1800,10 +1800,17 @@ void Neighbor::print_pairwise_info() else out += fmt::format(", half/full from ({})",rq->halffulllist+1); } else if (rq->skip) { - if (rq->trim) - out += fmt::format(", skip trim from ({})",rq->skiplist+1); - else - out += fmt::format(", skip from ({})",rq->skiplist+1); + if (rq->molskip) { + if (rq->trim) + out += fmt::format(", molskip trim from ({})",rq->skiplist+1); + else + out += fmt::format(", molskip from ({})",rq->skiplist+1); + } else { + if (rq->trim) + out += fmt::format(", skip trim from ({})",rq->skiplist+1); + else + out += fmt::format(", skip from ({})",rq->skiplist+1); + } } out += "\n"; diff --git a/src/neighbor.h b/src/neighbor.h index 8533fe5efa..5e0165b411 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -334,7 +334,8 @@ namespace NeighConst { NP_HALF_FULL = 1 << 23, NP_OFF2ON = 1 << 24, NP_MULTI_OLD = 1 << 25, - NP_TRIM = 1 << 26 + NP_TRIM = 1 << 26, + NP_INTRA = 1 << 27 }; enum { diff --git a/src/npair_halffull.h b/src/npair_halffull.h index 41d2e37dc4..99d734dba7 100644 --- a/src/npair_halffull.h +++ b/src/npair_halffull.h @@ -53,13 +53,13 @@ typedef NPairHalffull<1, 0, 0> NPairHalffullNewton; NPairStyle(halffull/newton/skip, NPairHalffullNewton, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | - NP_ORTHO | NP_SKIP); + NP_ORTHO | NP_SKIP | NP_INTRA); typedef NPairHalffull<1, 1, 0> NPairHalffullNewtonTri; NPairStyle(halffull/newton/skip/tri, NPairHalffullNewtonTri, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | - NP_TRI | NP_SKIP); + NP_TRI | NP_SKIP | NP_INTRA); typedef NPairHalffull<0, 0, 1> NPairHalffullTrimNewtoff; NPairStyle(halffull/trim/newtoff, @@ -71,7 +71,7 @@ typedef NPairHalffull<0, 0, 1> NPairHalffullTrimNewtoff; NPairStyle(halffull/trim/newtoff/skip, NPairHalffullTrimNewtoff, NP_HALF_FULL | NP_NEWTOFF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_HALF | - NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM); + NP_ORTHO | NP_TRI | NP_SKIP | NP_INTRA | NP_TRIM); typedef NPairHalffull<0, 0, 1> NPairHalffullTrimNewtoff; NPairStyle(halffull/trim/newtoff/ghost, @@ -83,7 +83,7 @@ typedef NPairHalffull<0, 0, 1> NPairHalffullTrimNewtoff; NPairStyle(halffull/trim/newtoff/skip/ghost, NPairHalffullTrimNewtoff, NP_HALF_FULL | NP_NEWTOFF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_HALF | - NP_ORTHO | NP_TRI | NP_SKIP | NP_GHOST | NP_TRIM); + NP_ORTHO | NP_TRI | NP_SKIP | NP_INTRA | NP_GHOST | NP_TRIM); typedef NPairHalffull<1, 0, 1> NPairHalffullTrimNewton; NPairStyle(halffull/trim/newton, @@ -101,13 +101,13 @@ typedef NPairHalffull<1, 0, 1> NPairHalffullTrimNewton; NPairStyle(halffull/trim/newton/skip, NPairHalffullTrimNewton, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | - NP_ORTHO | NP_SKIP | NP_TRIM); + NP_ORTHO | NP_SKIP | NP_INTRA | NP_TRIM); typedef NPairHalffull<1, 1, 1> NPairHalffullTrimNewtonTri; NPairStyle(halffull/trim/newton/tri/skip, NPairHalffullTrimNewtonTri, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | - NP_TRI | NP_SKIP | NP_TRIM); + NP_TRI | NP_SKIP | NP_INTRA | NP_TRIM); // clang-format on #else diff --git a/src/npair_skip.cpp b/src/npair_skip.cpp index 6afb43bc16..91d94edea0 100644 --- a/src/npair_skip.cpp +++ b/src/npair_skip.cpp @@ -17,6 +17,7 @@ #include "error.h" #include "my_page.h" #include "neigh_list.h" +#include "neigh_request.h" using namespace LAMMPS_NS; @@ -41,11 +42,13 @@ void NPairSkipTemp::build(NeighList *list) int *type = atom->type; int nlocal = atom->nlocal; + tagint *molecule = atom->molecule; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; + int molskip = list->molskip; int *ilist_skip = list->listskip->ilist; int *numneigh_skip = list->listskip->numneigh; @@ -71,7 +74,8 @@ void NPairSkipTemp::build(NeighList *list) for (ii = 0; ii < num_skip; ii++) { i = ilist_skip[ii]; itype = type[i]; - if (iskip[itype]) continue; + + if (!molskip && iskip[itype]) continue; if (TRIM) { xtmp = x[i][0]; @@ -90,8 +94,10 @@ void NPairSkipTemp::build(NeighList *list) for (jj = 0; jj < jnum; jj++) { joriginal = jlist[jj]; j = joriginal & NEIGHMASK; - if (ijskip[itype][type[j]]) continue; - + if (!molskip && ijskip[itype][type[j]]) continue; + if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue; + if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue; + if (TRIM) { delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/npair_skip.h b/src/npair_skip.h index cb0d201555..31c4910c93 100644 --- a/src/npair_skip.h +++ b/src/npair_skip.h @@ -16,41 +16,43 @@ typedef NPairSkipTemp<0> NPairSkip; NPairStyle(skip, NPairSkip, - NP_SKIP | NP_HALF | NP_FULL | + NP_SKIP | NP_HALF | NP_FULL | NP_INTRA | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI); typedef NPairSkipTemp<0> NPairSkip; NPairStyle(skip/ghost, NPairSkip, - NP_SKIP | NP_HALF | NP_FULL | + NP_SKIP | NP_HALF | NP_FULL | NP_INTRA | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_GHOST); typedef NPairSkipTemp<0> NPairSkipSize; NPairStyle(skip/half/size, NPairSkipSize, - NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | + NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_INTRA | + NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI); typedef NPairSkipTemp<1> NPairSkipTrim; NPairStyle(skip/trim, NPairSkipTrim, - NP_SKIP | NP_HALF | NP_FULL | + NP_SKIP | NP_HALF | NP_FULL | NP_INTRA | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM); typedef NPairSkipTemp<1> NPairSkipTrim; NPairStyle(skip/ghost/trim, NPairSkipTrim, - NP_SKIP | NP_HALF | NP_FULL | + NP_SKIP | NP_HALF | NP_FULL | NP_INTRA | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM); typedef NPairSkipTemp<1> NPairSkipTrimSize; NPairStyle(skip/trim/half/size, NPairSkipTrimSize, - NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | + NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_INTRA | + NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM); // clang-format on diff --git a/src/npair_skip_respa.h b/src/npair_skip_respa.h index af25e84faf..00d5b7f70b 100644 --- a/src/npair_skip_respa.h +++ b/src/npair_skip_respa.h @@ -16,14 +16,14 @@ typedef NPairSkipRespaTemp<0> NPairSkipRespa; NPairStyle(skip/half/respa, NPairSkipRespa, - NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | + NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | NP_INTRA | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI); typedef NPairSkipRespaTemp<1> NPairSkipTrimRespa; NPairStyle(skip/trim/half/respa, NPairSkipTrimRespa, - NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | + NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | NP_INTRA | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM); diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index d257973617..c6a1b4695f 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -321,11 +321,12 @@ void PairHybrid::settings(int narg, char **arg) iarg = 0; nstyles = 0; + const std::string mystyle = force->pair_style; while (iarg < narg) { if (utils::strmatch(arg[iarg],"^hybrid")) - error->all(FLERR,"Pair style hybrid cannot have hybrid as a sub-style"); + error->all(FLERR,"Pair style {} cannot have hybrid as a sub-style", mystyle); if (strcmp(arg[iarg],"none") == 0) - error->all(FLERR,"Pair style hybrid cannot have none as a sub-style"); + error->all(FLERR,"Pair style {} cannot have none as a sub-style", mystyle); styles[nstyles] = force->new_pair(arg[iarg],1,dummy); keywords[nstyles] = force->store_style(arg[iarg],0); @@ -345,6 +346,9 @@ void PairHybrid::settings(int narg, char **arg) nstyles++; } + if (utils::strmatch(mystyle,"^hybrid/molecular") && (nstyles != 2)) + error->all(FLERR, "Pair style {} must have exactly two sub-styles", mystyle); + delete[] cutmax_style; cutmax_style = new double[nstyles]; memset(cutmax_style, 0, nstyles*sizeof(double));