git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13757 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/fix_respa.cpp

@@ -12,6 +12,7 @@
 ------------------------------------------------------------------------- */
 
 #include "stdlib.h"
+#include "string.h"
 #include "fix_respa.h"
 #include "atom.h"
 #include "force.h"
@@ -29,10 +30,18 @@ FixRespa::FixRespa(LAMMPS *lmp, int narg, char **arg) :
 
   nlevels = force->inumeric(FLERR,arg[3]);
 
+  // optional arguments
+  store_torque = 0;
+  for (int iarg=4; iarg < narg; ++iarg) {
+    if (strcmp(arg[iarg],"torque") == 0)
+       store_torque = 1;
+  }
+
   // perform initial allocation of atom-based arrays
   // register with Atom class
 
   f_level = NULL;
+  t_level = NULL;
   grow_arrays(atom->nmax);
   atom->add_callback(0);
 }
@@ -48,6 +57,7 @@ FixRespa::~FixRespa()
   // delete locally stored arrays
 
   memory->destroy(f_level);
+  if (store_torque) memory->destroy(t_level);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -64,6 +74,7 @@ int FixRespa::setmask()
 double FixRespa::memory_usage()
 {
   double bytes = atom->nmax*nlevels*3 * sizeof(double);
+  if (store_torque) bytes += atom->nmax*nlevels*3 * sizeof(double);
   return bytes;
 }
 
@@ -74,6 +85,7 @@ double FixRespa::memory_usage()
 void FixRespa::grow_arrays(int nmax)
 {
   memory->grow(f_level,nmax,nlevels,3,"fix_respa:f_level");
+  if (store_torque) memory->grow(t_level,nmax,nlevels,3,"fix_respa:t_level");
 }
 
 /* ----------------------------------------------------------------------
@@ -87,6 +99,13 @@ void FixRespa::copy_arrays(int i, int j, int delflag)
     f_level[j][k][1] = f_level[i][k][1];
     f_level[j][k][2] = f_level[i][k][2];
   }
+  if (store_torque) {
+    for (int k = 0; k < nlevels; k++) {
+      t_level[j][k][0] = t_level[i][k][0];
+      t_level[j][k][1] = t_level[i][k][1];
+      t_level[j][k][2] = t_level[i][k][2];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -101,6 +120,13 @@ int FixRespa::pack_exchange(int i, double *buf)
     buf[m++] = f_level[i][k][1];
     buf[m++] = f_level[i][k][2];
   }
+  if (store_torque) {
+    for (int k = 0; k < nlevels; k++) {
+      buf[m++] = t_level[i][k][0];
+      buf[m++] = t_level[i][k][1];
+      buf[m++] = t_level[i][k][2];
+    }
+  }
   return m;
 }
 
@@ -116,5 +142,12 @@ int FixRespa::unpack_exchange(int nlocal, double *buf)
     f_level[nlocal][k][1] = buf[m++];
     f_level[nlocal][k][2] = buf[m++];
   }
+  if (store_torque) {
+    for (int k = 0; k < nlevels; k++) {
+      t_level[nlocal][k][0] = buf[m++];
+      t_level[nlocal][k][1] = buf[m++];
+      t_level[nlocal][k][2] = buf[m++];
+    }
+  }
   return m;
 }

src/fix_respa.h

@@ -43,7 +43,9 @@ class FixRespa : public Fix {
 
  private:
   int nlevels;
-  double ***f_level;            // force at each rRESPA level
+  int store_torque;  // 1 if torques should be stored in addition to forces
+  double ***f_level; // force at each rRESPA level
+  double ***t_level; // torque at each rRESPA level
 };
 
 }

src/pair_hybrid.cpp

@@ -25,6 +25,7 @@
 #include "comm.h"
 #include "memory.h"
 #include "error.h"
+#include "respa.h"
 
 using namespace LAMMPS_NS;
 
@@ -36,8 +37,11 @@ PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp)
   styles = NULL;
   keywords = NULL;
   multiple = NULL;
+  special_lj = NULL;
+  special_coul = NULL;
 
   outerflag = 0;
+  respaflag = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -45,13 +49,20 @@ PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp)
 PairHybrid::~PairHybrid()
 {
   if (nstyles) {
-    for (int m = 0; m < nstyles; m++) delete styles[m];
+    for (int m = 0; m < nstyles; m++) {
-    for (int m = 0; m < nstyles; m++) delete [] keywords[m];
+      delete styles[m];
+      delete [] keywords[m];
+      if (special_lj[m])   delete [] special_lj[m];
+      if (special_coul[m]) delete [] special_coul[m];
+    }
   }
   delete [] styles;
   delete [] keywords;
   delete [] multiple;
 
+  delete [] special_lj;
+  delete [] special_coul;
+
   delete [] svector;
 
   if (allocated) {
@@ -98,15 +109,36 @@ void PairHybrid::compute(int eflag, int vflag)
   if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4;
   else vflag_substyle = vflag;
 
+  double *saved_special = save_special();
+
+  // check if we are running with r-RESPA using the hybrid keyword
+
+  Respa *respa = NULL;
+  if (strstr(update->integrate_style,"respa")) {
+    respa = (Respa *) update->integrate;
+    if (respa->nhybrid_styles > 0) respaflag = 1;
+  }
+
   for (m = 0; m < nstyles; m++) {
 
-    // invoke compute() unless compute flag is turned off or
+    set_special(m);
-    // outerflag is set and sub-style has a compute_outer() method
 
-    if (styles[m]->compute_flag == 0) continue;
+    if (!respaflag || (respaflag && respa->hybrid_compute[m])) {
-    if (outerflag && styles[m]->respa_enable) 
+
-      styles[m]->compute_outer(eflag,vflag_substyle);
+      // invoke compute() unless compute flag is turned off or
-    else styles[m]->compute(eflag,vflag_substyle);
+      // outerflag is set and sub-style has a compute_outer() method
+
+      if (styles[m]->compute_flag == 0) continue;
+      if (outerflag && styles[m]->respa_enable) 
+        styles[m]->compute_outer(eflag,vflag_substyle);
+      else styles[m]->compute(eflag,vflag_substyle);
+    }
+
+    restore_special(saved_special);
+
+    // jump to next sub-style if r-RESPA does not want global accumulated data
+
+    if (respaflag && !respa->tally_global) continue;
 
     if (eflag_global) {
       eng_vdwl += styles[m]->eng_vdwl;
@@ -131,6 +163,8 @@ void PairHybrid::compute(int eflag, int vflag)
     }
   }
 
+  delete [] saved_special;
+
   if (vflag_fdotr) virial_fdotr_compute();
 }
 
@@ -215,6 +249,9 @@ void PairHybrid::settings(int narg, char **arg)
   keywords = new char*[narg];
   multiple = new int[narg];
 
+  special_lj = new double*[narg];
+  special_coul = new double*[narg];
+
   // allocate each sub-style
   // allocate uses suffix, but don't store suffix version in keywords,
   //   else syntax in coeff() will not match
@@ -233,6 +270,7 @@ void PairHybrid::settings(int narg, char **arg)
 
     styles[nstyles] = force->new_pair(arg[iarg],1,dummy);
     force->store_style(keywords[nstyles],arg[iarg],0);
+    special_lj[nstyles] = special_coul[nstyles] = NULL;
 
     jarg = iarg + 1;
     while (jarg < narg && !force->pair_map->count(arg[jarg])) jarg++;
@@ -421,6 +459,28 @@ void PairHybrid::init_style()
     if (used == 0) error->all(FLERR,"Pair hybrid sub-style is not used");
   }
 
+  // check if special_lj/special_coul overrides are compatible
+  for (istyle = 0; istyle < nstyles; istyle++) {
+    if (special_lj[istyle]) {
+      for (i = 1; i < 4; ++i) {
+        if (((force->special_lj[i] == 0.0) || (force->special_lj[i] == 1.0))
+            && (force->special_lj[i] != special_lj[istyle][i]))
+          error->all(FLERR,"Pair_modify special setting incompatible with"
+                     " global special_bonds setting");
+      }
+    }
+
+    if (special_coul[istyle]) {
+      for (i = 1; i < 4; ++i) {
+        if (((force->special_coul[i] == 0.0)
+             || (force->special_coul[i] == 1.0))
+            && (force->special_coul[i] != special_coul[istyle][i]))
+          error->all(FLERR,"Pair_modify special setting incompatible with"
+                     "global special_bonds setting");
+      }
+    }
+  }
+
   // each sub-style makes its neighbor list request(s)
 
   for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();
@@ -610,6 +670,13 @@ void PairHybrid::write_restart(FILE *fp)
     fwrite(&n,sizeof(int),1,fp);
     fwrite(keywords[m],sizeof(char),n,fp);
     styles[m]->write_restart_settings(fp);
+    // write out per style special settings, if present
+    n = (special_lj[m] == NULL) ? 0 : 1;
+    fwrite(&n,sizeof(int),1,fp);
+    if (n) fwrite(special_lj[m],sizeof(double),4,fp);
+    n = (special_coul[m] == NULL) ? 0 : 1;
+    fwrite(&n,sizeof(int),1,fp);
+    if (n) fwrite(special_coul[m],sizeof(double),4,fp);
   }
 }
 
@@ -629,6 +696,9 @@ void PairHybrid::read_restart(FILE *fp)
   keywords = new char*[nstyles];
   multiple = new int[nstyles];
 
+  special_lj = new double*[nstyles];
+  special_coul = new double*[nstyles];
+
   // each sub-style is created via new_pair()
   // each reads its settings, but no coeff info
 
@@ -641,6 +711,22 @@ void PairHybrid::read_restart(FILE *fp)
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
     styles[m] = force->new_pair(keywords[m],0,dummy);
     styles[m]->read_restart_settings(fp);
+    // read back per style special settings, if present
+    special_lj[m] = special_coul[m] = NULL;
+    if (me == 0) fread(&n,sizeof(int),1,fp);
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    if (n > 0 ) {
+      special_lj[m] = new double[4];
+      if (me == 0) fread(special_lj[m],sizeof(double),4,fp);
+      MPI_Bcast(special_lj[m],4,MPI_DOUBLE,0,world);
+    }
+    if (me == 0) fread(&n,sizeof(int),1,fp);
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    if (n > 0 ) {
+      special_coul[m] = new double[4];
+      if (me == 0) fread(special_coul[m],sizeof(double),4,fp);
+      MPI_Bcast(special_coul[m],4,MPI_DOUBLE,0,world);
+    }
   }
 
   // multiple[i] = 1 to M if sub-style used multiple times, else 0
@@ -681,6 +767,11 @@ double PairHybrid::single(int i, int j, int itype, int jtype,
       if (styles[map[itype][jtype][m]]->single_enable == 0)
         error->one(FLERR,"Pair hybrid sub-style does not support single call");
 
+      if ((special_lj[map[itype][jtype][m]] != NULL) ||
+          (special_coul[map[itype][jtype][m]] != NULL))
+        error->one(FLERR,"Pair hybrid single calls do not support"
+                   " per sub-style special bond values");
+
       esum += styles[map[itype][jtype][m]]->
         single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone);
       fforce += fone;
@@ -716,8 +807,14 @@ void PairHybrid::modify_params(int narg, char **arg)
       if (strcmp(arg[1],keywords[m]) == 0) break;
     if (m == nstyles) error->all(FLERR,"Unknown pair_modify hybrid sub-style");
     if (multiple[m] == 0) {
-      Pair::modify_params(narg-2,&arg[2]);
+      // augment special settings for this one pair style
-      styles[m]->modify_params(narg-2,&arg[2]);
+      if (strcmp(arg[2],"special") == 0) {
+        if (narg < 7) error->all(FLERR,"Illegal pair_modify special command");
+        modify_special(m,narg-3,&arg[3]);
+      } else {
+        Pair::modify_params(narg-2,&arg[2]);
+        styles[m]->modify_params(narg-2,&arg[2]);
+      }
     } else {
       if (narg < 3) error->all(FLERR,"Illegal pair_modify command");
       int multiflag = force->inumeric(FLERR,arg[2]);
@@ -725,13 +822,98 @@ void PairHybrid::modify_params(int narg, char **arg)
         if (strcmp(arg[1],keywords[m]) == 0 && multiflag == multiple[m]) break;
       if (m == nstyles) 
         error->all(FLERR,"Unknown pair_modify hybrid sub-style");
-      Pair::modify_params(narg-3,&arg[3]);
+      // augment special settings for this one pair style
-      styles[m]->modify_params(narg-3,&arg[3]);
+      if (strcmp(arg[3],"special") == 0) {
+        if (narg < 8) error->all(FLERR,"Illegal pair_modify special command");
+        modify_special(m,narg-4,&arg[4]);
+      } else {
+        Pair::modify_params(narg-3,&arg[3]);
+        styles[m]->modify_params(narg-3,&arg[3]);
+      }
     }
     
   } else {
-    Pair::modify_params(narg,arg);
+    // augment special settings for all pair styles
-    for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg);
+    if (strcmp(arg[0],"special") == 0) {
+      if (narg < 5) error->all(FLERR,"Illegal pair_modify special command");
+      for (int m = 0; m < nstyles; m++) modify_special(m,narg-1,&arg[1]);
+    } else {
+      Pair::modify_params(narg,arg);
+      for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg);
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   store a local per pair style override for special_lj and special_coul
+------------------------------------------------------------------------- */
+
+void PairHybrid::modify_special(int m, int narg, char **arg)
+{
+  double special[4];
+  int i;
+
+  special[0] = 1.0;
+  special[1] = force->numeric(FLERR,arg[1]);
+  special[2] = force->numeric(FLERR,arg[2]);
+  special[3] = force->numeric(FLERR,arg[3]);
+
+  if (strcmp(arg[0],"lj/coul") == 0) {
+    if (!special_lj[m]) special_lj[m] = new double[4];
+    if (!special_coul[m]) special_coul[m] = new double[4];
+    for (i = 0; i < 4; ++i)
+      special_lj[m][i] = special_coul[m][i] = special[i];
+
+  } else if (strcmp(arg[0],"lj") == 0) {
+    if (!special_lj[m]) special_lj[m] = new double[4];
+    for (i = 0; i < 4; ++i)
+      special_lj[m][i] = special[i];
+
+  } else if (strcmp(arg[0],"coul") == 0) {
+    if (!special_coul[m]) special_coul[m] = new double[4];
+    for (i = 0; i < 4; ++i)
+      special_coul[m][i] = special[i];
+
+  } else error->all(FLERR,"Illegal pair_modify special command");
+}
+
+/* ----------------------------------------------------------------------
+   override global special bonds settings with per substyle values
+------------------------------------------------------------------------- */
+
+void PairHybrid::set_special(int m)
+{
+  int i;
+  if (special_lj[m])
+    for (i = 0; i < 4; ++i) force->special_lj[i] = special_lj[m][i];
+  if (special_coul[m])
+    for (i = 0; i < 4; ++i) force->special_coul[i] = special_coul[m][i];
+}
+
+/* ----------------------------------------------------------------------
+   store global special settings
+------------------------------------------------------------------------- */
+
+double * PairHybrid::save_special()
+{
+  double *saved = new double[8];
+
+  for (int i = 0; i < 4; ++i) {
+    saved[i] = force->special_lj[i];
+    saved[i+4] = force->special_coul[i];
+  }
+  return saved;
+}
+
+/* ----------------------------------------------------------------------
+   restore global special settings from saved data
+------------------------------------------------------------------------- */
+
+void PairHybrid::restore_special(double *saved)
+{
+  for (int i = 0; i < 4; ++i) {
+    force->special_lj[i] = saved[i];
+    force->special_coul[i] = saved[i+4];
   }
 }
 

src/pair_hybrid.h

@@ -26,12 +26,11 @@ PairStyle(hybrid,PairHybrid)
 namespace LAMMPS_NS {
 
 class PairHybrid : public Pair {
+  friend class FixGPU;
+  friend class FixOMP;
+  friend class Force;
+  friend class Respa;
  public:
-  int nstyles;                  // # of sub-styles
-  Pair **styles;                // list of Pair style classes
-  char **keywords;              // style name of each Pair style
-  int *multiple;                // 0 if style used once, else Mth instance
-
   PairHybrid(class LAMMPS *);
   virtual ~PairHybrid();
   void compute(int, int);
@@ -54,13 +53,28 @@ class PairHybrid : public Pair {
   int check_ijtype(int, int, char *);
 
  protected:
+  int nstyles;                  // # of sub-styles
+  Pair **styles;                // list of Pair style classes
+  char **keywords;              // style name of each Pair style
+  int *multiple;                // 0 if style used once, else Mth instance
+
   int outerflag;                // toggle compute() when invoked by outer()
+  int respaflag;                // 1 if different substyles are assigned to
+                                // different r-RESPA levels
 
   int **nmap;                   // # of sub-styles itype,jtype points to
   int ***map;                   // list of sub-styles itype,jtype points to
+  double **special_lj;          // list of per style LJ exclusion factors
+  double **special_coul;        // list of per style Coulomb exclusion factors
 
   void allocate();
   void flags();
+
+  void modify_special(int, int, char**);
+  double *save_special();
+  void set_special(int);
+  void restore_special(double *);
+
   virtual void modify_requests();
 };
 

src/respa.cpp

@@ -38,6 +38,7 @@
 #include "timer.h"
 #include "memory.h"
 #include "error.h"
+#include "pair_hybrid.h"
 
 using namespace LAMMPS_NS;
 
@@ -65,6 +66,11 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
   level_pair = level_kspace = -1;
   level_inner = level_middle = level_outer = -1;
 
+  // defaults for hybrid pair styles
+  nhybrid_styles = 0;
+  tally_global = 1;
+  pair_compute = 1;
+
   int iarg = nlevels;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"bond") == 0) {
@@ -107,6 +113,22 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_kspace = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"hybrid") == 0) {
+      // the hybrid keyword requires a hybrid pair style
+      if (!strstr(force->pair_style,"hybrid"))
+        error->all(FLERR,"Illegal run_style respa command");
+      PairHybrid *hybrid = (PairHybrid *) force->pair;
+      nhybrid_styles = hybrid->nstyles;
+      // each hybrid sub-style needs to be assigned to a respa level
+      if (iarg+nhybrid_styles > narg)
+        error->all(FLERR,"Illegal run_style respa command");
+      hybrid_level = new int[nhybrid_styles];
+      hybrid_compute = new int[nhybrid_styles];
+      for (int i=0; i < nhybrid_styles; ++i) {
+        ++iarg;
+        hybrid_level[i] = force->inumeric(FLERR,arg[iarg])-1;
+      }
+      ++iarg;
     } else error->all(FLERR,"Illegal run_style respa command");
   }
 
@@ -127,6 +149,10 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
   if (level_middle >= 0 && level_inner == -1)
     error->all(FLERR,"Cannot set respa middle without inner/outer");
 
+  // cannot combine hybrid with any of pair/inner/middle/outer
+  if ((nhybrid_styles > 0) && (level_pair >= 0 || level_inner >= 0
+                               || level_middle >= 0 || level_outer >= 0))
+    error->all(FLERR,"Cannot set respa hybrid and any of pair/inner/middle/outer");
   // set defaults if user did not specify level
   // bond to innermost level
   // angle same as bond, dihedral same as angle, improper same as dihedral
@@ -138,9 +164,23 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
   if (level_angle == -1) level_angle = level_bond;
   if (level_dihedral == -1) level_dihedral = level_angle;
   if (level_improper == -1) level_improper = level_dihedral;
-  if (level_pair == -1 && level_inner == -1) level_pair = nlevels-1;
+
+  if (level_pair == -1 && level_inner == -1 && nhybrid_styles < 1)
+    level_pair = nlevels-1;
+
   if (level_kspace == -1 && level_pair >= 0) level_kspace = level_pair;
-  if (level_kspace == -1 && level_pair == -1) level_kspace = level_outer;
+  if (level_kspace == -1 && level_pair == -1) {
+    if (nhybrid_styles < 1) {
+      level_kspace = level_outer;
+    } else {
+      int max_hybrid_level = -1;
+      for (int i=0; i < nhybrid_styles; ++i) {
+        if (max_hybrid_level < hybrid_level[i])
+          max_hybrid_level = hybrid_level[i];
+      }
+      level_kspace = max_hybrid_level;
+    }
+  }
 
   // print respa levels
 
@@ -157,6 +197,9 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
         if (level_inner == i) fprintf(screen," pair-inner");
         if (level_middle == i) fprintf(screen," pair-middle");
         if (level_outer == i) fprintf(screen," pair-outer");
+        for (int j=0;j<nhybrid_styles;j++) {
+          if (hybrid_level[j] == i) fprintf(screen, " hybrid-%d",j+1);
+        }
         if (level_kspace == i) fprintf(screen," kspace");
         fprintf(screen,"\n");
       }
@@ -173,6 +216,9 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
         if (level_inner == i) fprintf(logfile," pair-inner");
         if (level_middle == i) fprintf(logfile," pair-middle");
         if (level_outer == i) fprintf(logfile," pair-outer");
+        for (int j=0;j<nhybrid_styles;j++) {
+          if (hybrid_level[j] == i) fprintf(logfile, " hybrid-%d",j+1);
+        }
         if (level_kspace == i) fprintf(logfile," kspace");
         fprintf(logfile,"\n");
       }
@@ -188,7 +234,7 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
     if (level_pair < level_improper || level_kspace < level_pair)
       error->all(FLERR,"Invalid order of forces within respa levels");
   }
-  if (level_pair == -1 && level_middle == -1) {
+  if (level_pair == -1 && level_middle == -1 && nhybrid_styles < 1) {
     if (level_inner < level_improper || level_outer < level_inner ||
         level_kspace < level_outer)
       error->all(FLERR,"Invalid order of forces within respa levels");
@@ -223,6 +269,12 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
     cutoff[3] = cutoff[1];
   }
 
+  // ensure that pair->compute() is run properly when the "hybrid" keyword is not used. 
+  if (nhybrid_styles < 1) {
+    pair_compute = 1;
+    tally_global = 1;
+  }
+
   // allocate other needed arrays
 
   newton = new int[nlevels];
@@ -236,6 +288,10 @@ Respa::~Respa()
   delete [] loop;
   delete [] newton;
   delete [] step;
+  if (nhybrid_styles > 0) {
+    delete [] hybrid_level;
+    delete [] hybrid_compute;
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -254,13 +310,19 @@ void Respa::init()
   // create fix needed for storing atom-based respa level forces
   // will delete it at end of run
 
-  char **fixarg = new char*[4];
+  char **fixarg = new char*[5];
   fixarg[0] = (char *) "RESPA";
   fixarg[1] = (char *) "all";
   fixarg[2] = (char *) "RESPA";
   fixarg[3] = new char[8];
   sprintf(fixarg[3],"%d",nlevels);
-  modify->add_fix(4,fixarg);
+  // if supported, we also store torques on a per-level basis
+  if (atom->torque_flag) {
+    fixarg[4] = (char *) "torque";
+    modify->add_fix(5,fixarg);
+  } else {
+    modify->add_fix(4,fixarg);
+  }
   delete [] fixarg[3];
   delete [] fixarg;
   fix_respa = (FixRespa *) modify->fix[modify->nfix-1];
@@ -309,6 +371,11 @@ void Respa::init()
       if (level_pair == ilevel || level_inner == ilevel ||
           level_middle == ilevel || level_outer == ilevel)
         newton[ilevel] = 1;
+
+      if (nhybrid_styles > 0) {
+        set_compute_flags(ilevel);
+        if (pair_compute) newton[ilevel] = 1;
+      }
     }
   }
 
@@ -360,6 +427,11 @@ void Respa::setup()
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     force_clear(newton[ilevel]);
     modify->setup_pre_force_respa(vflag,ilevel);
+
+    if (nhybrid_styles > 0) {
+      set_compute_flags(ilevel);
+      force->pair->compute(eflag,vflag);
+    }
     if (level_pair == ilevel && pair_compute_flag)
       force->pair->compute(eflag,vflag);
     if (level_inner == ilevel && pair_compute_flag)
@@ -428,6 +500,12 @@ void Respa::setup_minimal(int flag)
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     force_clear(newton[ilevel]);
     modify->setup_pre_force_respa(vflag,ilevel);
+
+    if (nhybrid_styles > 0) {
+      set_compute_flags(ilevel);
+      force->pair->compute(eflag,vflag);
+    }
+
     if (level_pair == ilevel && pair_compute_flag)
       force->pair->compute(eflag,vflag);
     if (level_inner == ilevel && pair_compute_flag)
@@ -574,6 +652,11 @@ void Respa::recurse(int ilevel)
       modify->pre_force_respa(vflag,ilevel,iloop);
 
     timer->stamp();
+    if (nhybrid_styles > 0) {
+      set_compute_flags(ilevel);
+      force->pair->compute(eflag,vflag);
+      timer->stamp(TIME_PAIR);
+    }
     if (level_pair == ilevel && pair_compute_flag) {
       force->pair->compute(eflag,vflag);
       timer->stamp(TIME_PAIR);
@@ -654,12 +737,19 @@ void Respa::copy_f_flevel(int ilevel)
 {
   double ***f_level = fix_respa->f_level;
   double **f = atom->f;
+  double ***t_level = fix_respa->t_level;
+  double **t = atom->torque;
   int n = atom->nlocal;
 
   for (int i = 0; i < n; i++) {
     f_level[i][ilevel][0] = f[i][0];
     f_level[i][ilevel][1] = f[i][1];
     f_level[i][ilevel][2] = f[i][2];
+    if (fix_respa->store_torque) {
+      t_level[i][ilevel][0] = t[i][0];
+      t_level[i][ilevel][1] = t[i][1];
+      t_level[i][ilevel][2] = t[i][2];
+    }
   }
 }
 
@@ -671,12 +761,19 @@ void Respa::copy_flevel_f(int ilevel)
 {
   double ***f_level = fix_respa->f_level;
   double **f = atom->f;
+  double ***t_level = fix_respa->t_level;
+  double **t = atom->torque;
   int n = atom->nlocal;
 
   for (int i = 0; i < n; i++) {
     f[i][0] = f_level[i][ilevel][0];
     f[i][1] = f_level[i][ilevel][1];
     f[i][2] = f_level[i][ilevel][2];
+    if (fix_respa->store_torque) {
+      t[i][0] = t_level[i][ilevel][0];
+      t[i][1] = t_level[i][ilevel][1];
+      t[i][2] = t_level[i][ilevel][2];
+    }
   }
 }
 
@@ -690,6 +787,8 @@ void Respa::sum_flevel_f()
 
   double ***f_level = fix_respa->f_level;
   double **f = atom->f;
+  double ***t_level = fix_respa->t_level;
+  double **t = atom->torque;
   int n = atom->nlocal;
 
   for (int ilevel = 1; ilevel < nlevels; ilevel++) {
@@ -697,6 +796,28 @@ void Respa::sum_flevel_f()
       f[i][0] += f_level[i][ilevel][0];
       f[i][1] += f_level[i][ilevel][1];
       f[i][2] += f_level[i][ilevel][2];
+      if (fix_respa->store_torque) {
+        t[i][0] += t_level[i][ilevel][0];
+        t[i][1] += t_level[i][ilevel][1];
+        t[i][2] += t_level[i][ilevel][2];
+      }
     }
   }
 }
+
+/*-----------------------------------------------------------------------
+  set flags for when some hybrid forces should be computed
+------------------------------------------------------------------------- */
+
+void Respa::set_compute_flags(int ilevel)
+{
+
+  if (nhybrid_styles < 1) return;
+
+  pair_compute = 0;
+  for (int i=0; i<nhybrid_styles; ++i) {
+    hybrid_compute[i] = (hybrid_level[i] == ilevel) ? 1 : 0;
+    if (hybrid_compute[i]) pair_compute = 1;
+  }
+  tally_global = (ilevel == nlevels-1) ? 1 : 0;
+}

src/respa.h

@@ -39,6 +39,12 @@ class Respa : public Integrate {
   int level_improper,level_pair,level_kspace;
   int level_inner,level_middle,level_outer;
 
+  int nhybrid_styles;     // number of hybrid pair styles
+  int *hybrid_level;      // level to compute pair hybrid sub-style at
+  int *hybrid_compute;    // selects whether to compute sub-style forces
+  int tally_global;       // 1 if pair style should tally global accumulators
+  int pair_compute;       // 1 if pair force need to be computed
+
   Respa(class LAMMPS *, int, char **);
   virtual ~Respa();
   virtual void init();
@@ -61,6 +67,7 @@ class Respa : public Integrate {
   virtual void recurse(int);
   void force_clear(int);
   void sum_flevel_f();
+  void set_compute_flags(int ilevel);
 };
 
 }
@@ -86,6 +93,12 @@ In the rRESPA integrator, you must compute pairwise potentials either
 all together (pair), or in pieces (inner/middle/outer).  You can't do
 both.
 
+E: Cannot set respa hybrid and any of pair/inner/middle/outer
+
+In the rRESPA integrator, you must compute pairwise potentials either
+all together (pair), with different cutoff regions (inner/middle/outer),
+or per hybrid sub-style (hybrid).  You cannot mix those.
+
 E: Must set both respa inner and outer
 
 Cannot use just the inner or outer option with respa without using the