initial version of pair style hybrid/molecular

doc/src/Commands_pair.rst

@@ -26,6 +26,7 @@ OPT.
    * :doc:`none <pair_none>`
    * :doc:`zero <pair_zero>`
    * :doc:`hybrid (k) <pair_hybrid>`
+   * :doc:`hybrid/molecular <pair_hybrid>`
    * :doc:`hybrid/overlay (k) <pair_hybrid>`
    * :doc:`hybrid/scaled <pair_hybrid>`
    * :doc:`kim <pair_kim>`
@@ -34,7 +35,6 @@ OPT.
    *
    *
    *
-   *
    * :doc:`adp (ko) <pair_adp>`
    * :doc:`agni (o) <pair_agni>`
    * :doc:`aip/water/2dm (t) <pair_aip_water_2dm>`

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
+The *hybrid*, *hybrid/overlay*, *hybrid/molecular*, and *hybrid/scaled*
-use of multiple pair styles in one simulation.  With the *hybrid* style,
+styles enable the use of multiple pair styles in one simulation.  With
-exactly one pair style is assigned to each pair of atom types.  With the
+the *hybrid* style, exactly one pair style is assigned to each pair of
-*hybrid/overlay* and *hybrid/scaled* styles, one or more pair styles can
+atom types.  With the *hybrid/overlay* and *hybrid/scaled* styles, one
-be assigned to each pair of atom types.  The assignment of pair styles
+or more pair styles can be assigned to each pair of atom types.  With
-to type pairs is made via the :doc:`pair_coeff <pair_coeff>` command.
+the hybrid/molecular style, pair styles are assigned to either intra-
-The major difference between the *hybrid/overlay* and *hybrid/scaled*
+or inter-molecular interactions.
-styles is that the *hybrid/scaled* adds a scale factor for each
+
-sub-style contribution to forces, energies and stresses.  Because of the
+The assignment of pair styles to type pairs is made via the
-added complexity, the *hybrid/scaled* style has more overhead and thus
+:doc:`pair_coeff <pair_coeff>` command.  The major difference between
-may be slower than *hybrid/overlay*.
+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.
 

doc/src/pair_style.rst

@@ -108,6 +108,7 @@ accelerated styles exist.
 
 * :doc:`none <pair_none>` - turn off pairwise interactions
 * :doc:`hybrid <pair_hybrid>` - multiple styles of pairwise interactions
+* :doc:`hybrid/molecular <pair_hybrid>` - different pair styles for intra- and inter-molecular interactions
 * :doc:`hybrid/overlay <pair_hybrid>` - multiple styles of superposed pairwise interactions
 * :doc:`hybrid/scaled <pair_hybrid>` - multiple styles of scaled superposed pairwise interactions
 * :doc:`zero <pair_zero>` - neighbor list but no interactions

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;
+    if (!molskip) {
-    iskip = new int[ntypes+1];
+      int ntypes = atom->ntypes;
-    memory->create(ijskip,ntypes+1,ntypes+1,"neigh_list:ijskip");
+      iskip = new int[ntypes+1];
-    int i,j;
+      memory->create(ijskip,ntypes+1,ntypes+1,"neigh_list:ijskip");
-    for (i = 1; i <= ntypes; i++) iskip[i] = nq->iskip[i];
+      int i,j;
-    for (i = 1; i <= ntypes; i++)
+      for (i = 1; i <= ntypes; i++) iskip[i] = nq->iskip[i];
-      for (j = 1; j <= ntypes; j++)
+      for (i = 1; i <= ntypes; i++)
-        ijskip[i][j] = nq->ijskip[i][j];
+        for (j = 1; j <= ntypes; j++)
+          ijskip[i][j] = nq->ijskip[i][j];
+    }
   }
 
   if (nq->halffull)

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
 

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;

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();

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)
+      if (rq->molskip) {
-        out += fmt::format(", skip trim from ({})",rq->skiplist+1);
+        if (rq->trim)
-      else
+          out += fmt::format(", molskip trim from ({})",rq->skiplist+1);
-        out += fmt::format(", skip 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";
 

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 {

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
 

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<TRIM>::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<int> *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<TRIM>::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<TRIM>::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];

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

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);
 

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));