Merge pull request #1150 from mkanski/compute_pair_multiple_pstyles

Extend compute pair to handle multiple instances of a sub-style in pair style hybrid

doc/src/compute_pair.txt

@@ -10,17 +10,20 @@ compute pair command :h3
 
 [Syntax:]
 
-compute ID group-ID pair pstyle evalue :pre
+compute ID group-ID pair pstyle \[nstyle\] \[evalue\]  :pre
 
-ID, group-ID are documented in "compute"_compute.html command
-pair = style name of this compute command
-pstyle = style name of a pair style that calculates additional values
-evalue = {epair} or {evdwl} or {ecoul} or blank (optional setting) :ul
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+pair = style name of this compute command :l
+pstyle = style name of a pair style that calculates additional values :l
+nsub = {n}-instance of a substyle, if a pair style is used multiple times in a hybrid style :l
+{evalue} = {epair} or {evdwl} or {ecoul} or blank (optional) :l
+:ule
 
 [Examples:]
 
 compute 1 all pair gauss
 compute 1 all pair lj/cut/coul/cut ecoul
+compute 1 all pair tersoff 2 epair
 compute 1 all pair reax :pre
 
 [Description:]
@@ -33,15 +36,19 @@ NOTE: The group specified for this command is [ignored].
 
 The specified {pstyle} must be a pair style used in your simulation
 either by itself or as a sub-style in a "pair_style hybrid or
-hybrid/overlay"_pair_hybrid.html command.
+hybrid/overlay"_pair_hybrid.html command. If the sub-style is
+used more than once, an additional number {nsub} has to be specified
+in order to choose which instance of the sub-style will be used by
+the compute. Not specifying the number in this case will cause the
+compute to fail.
 
-The {evalue} setting is optional; it may be left off the command.  All
+The {evalue} setting is optional.  All
 pair styles tally a potential energy {epair} which may be broken into
 two parts: {evdwl} and {ecoul} such that {epair} = {evdwl} + {ecoul}.
 If the pair style calculates Coulombic interactions, their energy will
 be tallied in {ecoul}.  Everything else (whether it is a Lennard-Jones
 style van der Waals interaction or not) is tallied in {evdwl}.  If
-{evalue} is specified as {epair} or left out, then {epair} is stored
+{evalue} is blank or specified as {epair}, then {epair} is stored
 as a global scalar by this compute.  This is useful when using
 "pair_style hybrid"_pair_hybrid.html if you want to know the portion
 of the total energy contributed by one sub-style.  If {evalue} is
@@ -82,4 +89,4 @@ the doc page for the pair style for details.
 
 [Default:]
 
-The default for {evalue} is {epair}.
+The keyword defaults are {evalue} = {epair}, nsub = 0.

src/compute_pair.cpp

@@ -13,6 +13,7 @@
 
 #include <mpi.h>
 #include <cstring>
+#include <cctype>
 #include "compute_pair.h"
 #include "update.h"
 #include "force.h"
@@ -29,7 +30,7 @@ ComputePair::ComputePair(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
   pstyle(NULL), pair(NULL), one(NULL)
 {
-  if (narg < 4 || narg > 5) error->all(FLERR,"Illegal compute pair command");
+  if (narg < 4) error->all(FLERR,"Illegal compute pair command");
 
   scalar_flag = 1;
   extscalar = 1;
@@ -41,19 +42,34 @@ ComputePair::ComputePair(LAMMPS *lmp, int narg, char **arg) :
   pstyle = new char[n];
   strcpy(pstyle,arg[3]);
 
-  if (narg == 5) {
-    if (strcmp(arg[4],"epair") == 0) evalue = EPAIR;
-    if (strcmp(arg[4],"evdwl") == 0) evalue = EVDWL;
-    if (strcmp(arg[4],"ecoul") == 0) evalue = ECOUL;
-  } else evalue = EPAIR;
+  int iarg = 4;
+  nsub = 0;
+  evalue = EPAIR;
+
+  if (narg > iarg) {
+    if (isdigit(arg[iarg][0])) {
+      nsub = force->inumeric(FLERR,arg[iarg]);
+      ++iarg;
+      if (nsub <= 0)
+        error->all(FLERR,"Illegal compute pair command");
+    }
+  }
+
+  if  (narg > iarg) {
+    if (strcmp(arg[iarg],"epair") == 0) evalue = EPAIR;
+    else if (strcmp(arg[iarg],"evdwl") == 0) evalue = EVDWL;
+    else if (strcmp(arg[iarg],"ecoul") == 0) evalue = ECOUL;
+    else error->all(FLERR, "Illegal compute pair command");
+    ++iarg;
+  }
 
   // check if pair style with and without suffix exists
 
-  pair = force->pair_match(pstyle,1);
+  pair = force->pair_match(pstyle,1,nsub);
   if (!pair && lmp->suffix) {
     strcat(pstyle,"/");
     strcat(pstyle,lmp->suffix);
-    pair = force->pair_match(pstyle,1);
+    pair = force->pair_match(pstyle,1,nsub);
   }
 
   if (!pair)
@@ -84,7 +100,7 @@ void ComputePair::init()
 {
   // recheck for pair style in case it has been deleted
 
-  pair = force->pair_match(pstyle,1);
+  pair = force->pair_match(pstyle,1,nsub);
   if (!pair)
     error->all(FLERR,"Unrecognized pair style in compute pair command");
 }

src/compute_pair.h

@@ -33,7 +33,7 @@ class ComputePair : public Compute {
   void compute_vector();
 
  private:
-  int evalue,npair;
+  int evalue,npair,nsub;
   char *pstyle;
   class Pair *pair;
   double *one;