Merge pull request #1511 from jdevemy/master

Add keyword hybrid/pair for compute_pressure

doc/src/compute_pressure.txt

@@ -16,12 +16,13 @@ ID, group-ID are documented in "compute"_compute.html command
 pressure = style name of this compute command
 temp-ID = ID of compute that calculates temperature, can be NULL if not needed
 zero or more keywords may be appended
-keyword = {ke} or {pair} or {bond} or {angle} or {dihedral} or {improper} or {kspace} or {fix} or {virial} :ul
+keyword = {ke} or {pair} or {bond} or {angle} or {dihedral} or {improper} or {kspace} or {fix} or {virial} or {pair/hybrid} :ul
 
 [Examples:]
 
 compute 1 all pressure thermo_temp
-compute 1 all pressure NULL pair bond :pre
+compute 1 all pressure NULL pair bond
+compute 1 all pressure NULL pair/hybrid lj/cut :pre
 
 [Description:]
 
@@ -67,6 +68,9 @@ extra keywords are listed, then only those components are summed to
 compute temperature or ke and/or the virial.  The {virial} keyword
 means include all terms except the kinetic energy {ke}.
 
+The {pair/hybrid} keyword means to only include contribution
+from a sub-style in a {hybrid} or {hybrid/overlay} pair style.
+
 Details of how LAMMPS computes the virial efficiently for the entire
 system, including for many-body potentials and accounting for the
 effects of periodic boundary conditions are discussed in

src/compute_pressure.cpp

@@ -21,6 +21,7 @@
 #include "fix.h"
 #include "force.h"
 #include "pair.h"
+#include "pair_hybrid.h"
 #include "bond.h"
 #include "angle.h"
 #include "dihedral.h"
@@ -65,6 +66,7 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
 
   // process optional args
 
+  pairhybridflag = 0;
   if (narg == 4) {
     keflag = 1;
     pairflag = 1;
@@ -78,6 +80,37 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
     int iarg = 4;
     while (iarg < narg) {
       if (strcmp(arg[iarg],"ke") == 0) keflag = 1;
+      else if (strcmp(arg[iarg],"pair/hybrid") == 0) {
+        int n = strlen(arg[++iarg]) + 1;
+        if (lmp->suffix) n += strlen(lmp->suffix) + 1;
+        pstyle = new char[n];
+        strcpy(pstyle,arg[iarg++]);
+
+        nsub = 0;
+
+        if (narg > iarg) {
+          if (isdigit(arg[iarg][0])) {
+            nsub = force->inumeric(FLERR,arg[iarg]);
+            ++iarg;
+            if (nsub <= 0)
+              error->all(FLERR,"Illegal compute pressure command");
+          }
+        }
+
+        // check if pair style with and without suffix exists
+
+        pairhybrid = (Pair *) force->pair_match(pstyle,1,nsub);
+        if (!pairhybrid && lmp->suffix) {
+          strcat(pstyle,"/");
+          strcat(pstyle,lmp->suffix);
+          pairhybrid = (Pair *) force->pair_match(pstyle,1,nsub);
+        }
+
+        if (!pairhybrid)
+          error->all(FLERR,"Unrecognized pair style in compute pressure command");
+
+        pairhybridflag = 1;
+      }
       else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
       else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
       else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
@@ -132,6 +165,20 @@ void ComputePressure::init()
     temperature = modify->compute[icompute];
   }
 
+  // recheck if pair style with and without suffix exists
+
+  if (pairhybridflag) {
+    pairhybrid = (Pair *) force->pair_match(pstyle,1,nsub);
+    if (!pairhybrid && lmp->suffix) {
+      strcat(pstyle,"/");
+      strcat(pstyle,lmp->suffix);
+      pairhybrid = (Pair *) force->pair_match(pstyle,1,nsub);
+    }
+
+    if (!pairhybrid)
+      error->all(FLERR,"Unrecognized pair style in compute pressure command");
+  }
+
   // detect contributions to virial
   // vptr points to all virial[6] contributions
 
@@ -139,6 +186,7 @@ void ComputePressure::init()
   nvirial = 0;
   vptr = NULL;
 
+  if (pairhybridflag && force->pair) nvirial++;
   if (pairflag && force->pair) nvirial++;
   if (bondflag && atom->molecular && force->bond) nvirial++;
   if (angleflag && atom->molecular && force->angle) nvirial++;
@@ -151,6 +199,11 @@ void ComputePressure::init()
   if (nvirial) {
     vptr = new double*[nvirial];
     nvirial = 0;
+    if (pairhybridflag && force->pair) {
+      PairHybrid *ph = (PairHybrid *) force->pair;
+      ph->no_virial_fdotr_compute = 1;
+      vptr[nvirial++] = pairhybrid->virial;
+    }
     if (pairflag && force->pair) vptr[nvirial++] = force->pair->virial;
     if (bondflag && force->bond) vptr[nvirial++] = force->bond->virial;
     if (angleflag && force->angle) vptr[nvirial++] = force->angle->virial;

src/compute_pressure.h

@@ -41,10 +41,16 @@ class ComputePressure : public Compute {
   Compute *temperature;
   char *id_temp;
   double virial[6];
+  int pairhybridflag;
+  class Pair *pairhybrid;
   int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag;
   int fixflag,kspaceflag;
 
   void virial_compute(int, int);
+
+ private:
+  char *pstyle;
+  int nsub;
 };
 
 }