Adding files for BPM

src/.gitignore

@@ -251,6 +251,23 @@
 /pair_mesont_tpm.cpp
 /pair_mesont_tpm.h
 
+/atom_vec_sphere_bpm.cpp
+/atom_vec_sphere_bpm.h
+/bond_bpm.cpp
+/bond_bpm.h
+/bond_bpm_rotational.cpp
+/bond_bpm_rotational.h
+/bond_bpm_simple.cpp
+/bond_bpm_simple.h
+/compute_nbond_atom.cpp
+/compute_nbond_atom.h
+/fix_bond_history.cpp
+/fix_bond_history.h
+/fix_nve_sphere_bpm.cpp
+/fix_nve_sphere_bpm.h
+/fix_update_special_bonds.cpp
+/fix_update_special_bonds.h
+
 /compute_adf.cpp
 /compute_adf.h
 /compute_contact_atom.cpp
@@ -779,8 +796,6 @@
 /fix_orient_eco.h
 /fix_orient_fcc.cpp
 /fix_orient_fcc.h
-/fix_pair_tracker.cpp
-/fix_pair_tracker.h
 /fix_peri_neigh.cpp
 /fix_peri_neigh.h
 /fix_phonon.cpp

src/BPM/atom_vec_sphere_bpm.cpp

@@ -0,0 +1,254 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "atom_vec_sphere_bpm.h"
+
+#include "atom.h"
+#include "error.h"
+#include "fix.h"
+#include "fix_adapt.h"
+#include "math_const.h"
+#include "modify.h"
+#include "utils.h"
+
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecSphereBPM::AtomVecSphereBPM(LAMMPS *lmp) : AtomVec(lmp)
+{
+  mass_type = PER_ATOM;
+  molecular = Atom::MOLECULAR;
+  bonds_allow = 1;
+
+  atom->molecule_flag = 1;
+  atom->sphere_flag = 1;
+  atom->radius_flag = atom->rmass_flag = atom->omega_flag =
+    atom->torque_flag = atom->quat_flag = 1;
+
+  // strings with peratom variables to include in each AtomVec method
+  // strings cannot contain fields in corresponding AtomVec default strings
+  // order of fields in a string does not matter
+  // except: fields_data_atom & fields_data_vel must match data file
+
+  fields_grow = (char *)
+    "molecule num_bond bond_type bond_atom nspecial special radius rmass omega torque quat";
+  fields_copy = (char *)
+    "molecule num_bond bond_type bond_atom nspecial special radius rmass omega quat";
+  fields_comm = (char *) "";
+  fields_comm_vel = (char *) "omega quat";
+  fields_reverse = (char *) "torque";
+  fields_border = (char *) "molecule radius rmass";
+  fields_border_vel = (char *) "molecule radius rmass omega quat";
+  fields_exchange = (char *)
+    "molecule num_bond bond_type bond_atom nspecial special radius rmass omega quat";
+  fields_restart = (char *)
+    "molecule num_bond bond_type bond_atom radius rmass omega quat";
+  fields_create = (char *)
+    "molecule num_bond nspecial radius rmass omega quat";
+  fields_data_atom = (char *) "id molecule type radius rmass x";
+  fields_data_vel = (char *) "id v omega";
+
+  bond_per_atom = 0;
+  bond_negative = NULL;
+}
+
+/* ----------------------------------------------------------------------
+   process sub-style args
+   optional arg = 0/1 for static/dynamic particle radii
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::process_args(int narg, char **arg)
+{
+  if (narg != 0 && narg != 1)
+    error->all(FLERR,"Illegal atom_style sphere/bpm command");
+
+  radvary = 0;
+  if (narg == 1) {
+    radvary = utils::numeric(FLERR,arg[0],true,lmp);
+    if (radvary < 0 || radvary > 1)
+      error->all(FLERR,"Illegal atom_style sphere/bpm command");
+  }
+
+  // dynamic particle radius and mass must be communicated every step
+
+  if (radvary) {
+    fields_comm = (char *) "radius rmass";
+    fields_comm_vel = (char *) "radius rmass omega";
+  }
+
+  // delay setting up of fields until now
+
+  setup_fields();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::init()
+{
+  AtomVec::init();
+
+  // check if optional radvary setting should have been set to 1
+
+  for (int i = 0; i < modify->nfix; i++)
+    if (strcmp(modify->fix[i]->style,"adapt") == 0) {
+      FixAdapt *fix = (FixAdapt *) modify->fix[i];
+      if (fix->diamflag && radvary == 0)
+        error->all(FLERR,"Fix adapt changes particle radii "
+                   "but atom_style sphere is not dynamic");
+    }
+}
+
+/* ----------------------------------------------------------------------
+   set local copies of all grow ptrs used by this class, except defaults
+   needed in replicate when 2 atom classes exist and it calls pack_restart()
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::grow_pointers()
+{
+  radius = atom->radius;
+  rmass = atom->rmass;
+  omega = atom->omega;
+  quat = atom->quat;
+
+  num_bond = atom->num_bond;
+  bond_type = atom->bond_type;
+  nspecial = atom->nspecial;
+}
+
+
+/* ----------------------------------------------------------------------
+   initialize non-zero atom quantities
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::create_atom_post(int ilocal)
+{
+  radius[ilocal] = 0.5;
+  rmass[ilocal] = 4.0*MY_PI/3.0 * 0.5*0.5*0.5;
+
+  quat[ilocal][0] = 1.0;
+  quat[ilocal][1] = 0.0;
+  quat[ilocal][2] = 0.0;
+  quat[ilocal][3] = 0.0;
+
+}
+
+/* ----------------------------------------------------------------------
+   modify values for AtomVec::pack_restart() to pack
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::pack_restart_pre(int ilocal)
+{
+  // insure bond_negative vector is needed length
+
+  if (bond_per_atom < atom->bond_per_atom) {
+    delete [] bond_negative;
+    bond_per_atom = atom->bond_per_atom;
+    bond_negative = new int[bond_per_atom];
+  }
+
+  // flip any negative types to positive and flag which ones
+
+  any_bond_negative = 0;
+  for (int m = 0; m < num_bond[ilocal]; m++) {
+    if (bond_type[ilocal][m] < 0) {
+      bond_negative[m] = 1;
+      bond_type[ilocal][m] = -bond_type[ilocal][m];
+      any_bond_negative = 1;
+    } else bond_negative[m] = 0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   unmodify values packed by AtomVec::pack_restart()
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::pack_restart_post(int ilocal)
+{
+  // restore the flagged types to their negative values
+
+  if (any_bond_negative) {
+    for (int m = 0; m < num_bond[ilocal]; m++)
+      if (bond_negative[m]) bond_type[ilocal][m] = -bond_type[ilocal][m];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   initialize other atom quantities after AtomVec::unpack_restart()
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::unpack_restart_init(int ilocal)
+{
+  nspecial[ilocal][0] = 0;
+  nspecial[ilocal][1] = 0;
+  nspecial[ilocal][2] = 0;
+}
+
+/* ----------------------------------------------------------------------
+   modify what AtomVec::data_atom() just unpacked
+   or initialize other atom quantities
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::data_atom_post(int ilocal)
+{
+  radius_one = 0.5 * atom->radius[ilocal];
+  radius[ilocal] = radius_one;
+  if (radius_one > 0.0)
+    rmass[ilocal] *= 4.0*MY_PI/3.0 * radius_one*radius_one*radius_one;
+
+  if (rmass[ilocal] <= 0.0)
+    error->one(FLERR,"Invalid density in Atoms section of data file");
+
+  omega[ilocal][0] = 0.0;
+  omega[ilocal][1] = 0.0;
+  omega[ilocal][2] = 0.0;
+
+  quat[ilocal][0] = 1.0;
+  quat[ilocal][1] = 0.0;
+  quat[ilocal][2] = 0.0;
+  quat[ilocal][3] = 0.0;
+
+  num_bond[ilocal] = 0;
+  nspecial[ilocal][0] = 0;
+  nspecial[ilocal][1] = 0;
+  nspecial[ilocal][2] = 0;
+}
+
+/* ----------------------------------------------------------------------
+   modify values for AtomVec::pack_data() to pack
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::pack_data_pre(int ilocal)
+{
+  radius_one = radius[ilocal];
+  rmass_one = rmass[ilocal];
+
+  radius[ilocal] *= 2.0;
+  if (radius_one!= 0.0)
+    rmass[ilocal] =
+      rmass_one / (4.0*MY_PI/3.0 * radius_one*radius_one*radius_one);
+}
+
+/* ----------------------------------------------------------------------
+   unmodify values packed by AtomVec::pack_data()
+------------------------------------------------------------------------- */
+
+void AtomVecSphereBPM::pack_data_post(int ilocal)
+{
+  radius[ilocal] = radius_one;
+  rmass[ilocal] = rmass_one;
+}

src/BPM/atom_vec_sphere_bpm.h

@@ -0,0 +1,83 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef ATOM_CLASS
+// clang-format off
+AtomStyle(sphere/bpm,AtomVecSphereBPM)
+// clang-format on
+#else
+
+#ifndef LMP_ATOM_VEC_SPHERE_BPM_H
+#define LMP_ATOM_VEC_SPHERE_BPM_H
+
+#include "atom_vec.h"
+
+namespace LAMMPS_NS {
+
+class AtomVecSphereBPM : public AtomVec {
+ public:
+  AtomVecSphereBPM(class LAMMPS *);
+  void process_args(int, char **);
+  void init();
+
+  void grow_pointers();
+  void create_atom_post(int);
+  void pack_restart_pre(int);
+  void pack_restart_post(int);
+  void unpack_restart_init(int);
+  void data_atom_post(int);
+  void pack_data_pre(int);
+  void pack_data_post(int);
+
+
+ private:
+  int *num_bond;
+  int **bond_type;
+  int **nspecial;
+
+  double *radius,*rmass;
+  double **omega, **torque, **quat;
+
+  int any_bond_negative;
+  int bond_per_atom;
+  int *bond_negative;
+
+  int radvary;
+  double radius_one,rmass_one;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Per-processor system is too big
+
+The number of owned atoms plus ghost atoms on a single
+processor must fit in 32-bit integer.
+
+E: Invalid atom type in Atoms section of data file
+
+Atom types must range from 1 to specified # of types.
+
+E: Invalid radius in Atoms section of data file
+
+Radius must be >= 0.0.
+
+E: Invalid density in Atoms section of data file
+
+Density value cannot be <= 0.0.
+
+*/

src/BPM/bond_bpm.cpp

@@ -0,0 +1,349 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "bond_bpm.h"
+
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "fix_store_local.h"
+#include "fix_update_special_bonds.h"
+#include "force.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp)
+{
+  id_fix_store_local = nullptr;
+  id_fix_prop_atom = nullptr;
+  fix_store_local = nullptr;
+  fix_update_special_bonds = nullptr;
+  
+  prop_atom_flag = 0;  
+  nvalues = 0;
+  output_data = nullptr;  
+  pack_choice = nullptr;
+  
+  r0_max_estimate = 0.0;
+  max_stretch = 1.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+BondBPM::~BondBPM()
+{
+  delete [] pack_choice;
+  delete [] id_fix_store_local;
+  delete [] id_fix_prop_atom;
+  memory->destroy(output_data);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::init_style()
+{
+  int ifix;
+  if (id_fix_store_local) {
+    ifix = modify->find_fix(id_fix_store_local);
+    if (ifix < 0) error->all(FLERR, "Cannot find fix store/local");
+    if (strcmp(modify->fix[ifix]->style, "store/local") != 0)
+      error->all(FLERR, "Incorrect fix style matched, not store/local");    
+    fix_store_local = (FixStoreLocal *) modify->fix[ifix];
+    fix_store_local->nvalues = nvalues;
+  } 
+  
+  ifix = modify->find_fix_by_style("update/special/bonds");
+  if (ifix >= 0)
+    fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix];
+  else 
+    fix_update_special_bonds = nullptr;
+
+
+  if (force->angle || force->dihedral || force->improper)
+    error->all(FLERR,
+               "Bond style bpm cannot be used with 3,4-body interactions");
+  if (atom->molecular == 2)
+    error->all(FLERR,
+               "Bond style bpm cannot be used with atom style template");
+               
+  // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists
+  if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 ||
+    force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0)
+    error->all(FLERR,"Bond style bpm requires 1-3 and 1-4 special weights of 1.0");    
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+   All args before store/local command are saved for potential args
+     for specific bond BPM substyles
+   All args after optional store/local command are variables stored
+     in the compute store/local
+------------------------------------------------------------------------- */
+
+void BondBPM::settings(int narg, char **arg)
+{ 
+  int iarg = 0;
+  while (iarg < narg) {
+    if (id_fix_store_local) {
+      if (strcmp(arg[iarg], "id1") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_id1;
+      } else if (strcmp(arg[iarg], "id2") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_id2;
+      } else if (strcmp(arg[iarg], "time") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_time;
+      } else if (strcmp(arg[iarg], "x") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_x;
+      } else if (strcmp(arg[iarg], "y") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_y;
+      } else if (strcmp(arg[iarg], "z") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_z;
+      } else if (strcmp(arg[iarg], "x/ref") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_x_ref;
+        prop_atom_flag = 1;      
+      } else if (strcmp(arg[iarg], "y/ref") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_y_ref;
+        prop_atom_flag = 1;      
+      } else if (strcmp(arg[iarg], "z/ref") == 0) {
+        pack_choice[nvalues++] = &BondBPM::pack_z_ref;
+        prop_atom_flag = 1;
+      } else {
+        error->all(FLERR, "Illegal bond_style command");
+      }
+    } else if (strcmp(arg[iarg], "store/local") == 0) {
+      id_fix_store_local = utils::strdup(arg[iarg+1]);
+      iarg ++;
+      nvalues = 0;      
+      pack_choice = new FnPtrPack[narg - iarg - 1];
+    }
+    iarg ++;
+  }
+  
+  if (id_fix_store_local) {
+    if (nvalues == 0) error->all(FLERR, 
+        "Bond style bpm/rotational must include at least one value to output");  
+    memory->create(output_data, nvalues, "bond/bpm:output_data");
+    
+    // Use store property to save reference positions as it can transfer to ghost atoms
+    if (prop_atom_flag == 1) {
+      
+      id_fix_prop_atom = utils::strdup("BPM_PROPERTY_ATOM");
+      int ifix = modify->find_fix(id_fix_prop_atom);
+      if (ifix < 0) {
+        modify->add_fix(fmt::format("{} all property/atom "
+            "d_BPM_X_REF d_BPM_Y_REF d_BPM_Z_REF ghost yes", id_fix_prop_atom));
+        ifix = modify->find_fix(id_fix_prop_atom);
+      }
+      
+      int type_flag;
+      int col_flag;
+      index_x_ref = atom->find_custom("BPM_X_REF", type_flag, col_flag);
+      index_y_ref = atom->find_custom("BPM_Y_REF", type_flag, col_flag);
+      index_z_ref = atom->find_custom("BPM_Z_REF", type_flag, col_flag);
+    
+      if (modify->fix[ifix]->restart_reset) {
+          modify->fix[ifix]->restart_reset = 0;
+      } else {
+        double *x_ref = atom->dvector[index_x_ref];
+        double *y_ref = atom->dvector[index_y_ref];
+        double *z_ref = atom->dvector[index_z_ref];
+    
+        double **x = atom->x;
+        for (int i = 0; i < atom->nlocal; i++) {
+          x_ref[i] = x[i][0];
+          y_ref[i] = x[i][1];
+          z_ref[i] = x[i][2];
+        }
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   used to check bond communiction cutoff - not perfect, estimates based on local-local only
+------------------------------------------------------------------------- */
+
+double BondBPM::equilibrium_distance(int i)
+{
+  // Ghost atoms may not yet be communicated, this may only be an estimate
+  if (r0_max_estimate == 0) {
+    int type, j;
+    double delx, dely, delz, r;
+    double **x = atom->x;
+    for (int i = 0; i < atom->nlocal; i ++) {
+      for (int m = 0; m < atom->num_bond[i]; m ++) {
+        type = atom->bond_type[i][m];
+        if (type == 0) continue;
+
+        j = atom->map(atom->bond_atom[i][m]);
+        if(j == -1) continue;
+
+        delx = x[i][0] - x[j][0];
+        dely = x[i][1] - x[j][1];
+        delz = x[i][2] - x[j][2];
+        domain->minimum_image(delx, dely, delz);
+        
+        r = sqrt(delx*delx + dely*dely + delz*delz);
+        if(r > r0_max_estimate) r0_max_estimate = r;
+      }
+    }
+
+    double temp;
+    MPI_Allreduce(&r0_max_estimate,&temp,1,MPI_DOUBLE,MPI_MAX,world);
+    r0_max_estimate = temp;
+
+    //if (comm->me == 0)
+    //  utils::logmesg(lmp,fmt::format("Estimating longest bond = {}\n",r0_max_estimate));
+  }
+
+  // Divide out heuristic prefactor added in comm class
+  return max_stretch*r0_max_estimate/1.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::process_broken(int i, int j)
+{
+  if (fix_store_local) {
+    for (int n = 0; n < nvalues; n++)
+      (this->*pack_choice[n])(n, i, j);
+
+    fix_store_local->add_data(output_data, i, j);
+  }
+  
+  if (fix_update_special_bonds) 
+    fix_update_special_bonds->add_broken_bond(i, j);
+
+  // Manually search and remove from atom arrays
+  // need to remove in case special bonds arrays rebuilt
+  int m, n;
+  int nlocal = atom->nlocal;
+  
+  tagint *tag = atom->tag;
+  tagint **bond_atom = atom->bond_atom;
+  int **bond_type = atom->bond_type;
+  int *num_bond = atom->num_bond;
+
+  if (i < nlocal) {
+    for (m = 0; m < num_bond[i]; m++) {
+   //   if(tag[i] == 25 or tag[j] == 25) printf("\t 1 ATOM LOOP %d-%d, %d/%d \n", tag[i], bond_atom[i][m], m, num_bond[i]);
+      if (bond_atom[i][m] == tag[j]) {
+  //      if(tag[i] == 25 or tag[j] == 25) printf("\t DELETED\n");
+        bond_type[i][m] = 0;
+        n = num_bond[i];
+        bond_type[i][m] = bond_type[i][n-1];
+        bond_atom[i][m] = bond_atom[i][n-1];
+        num_bond[i]--;
+        break;
+      }
+    }
+  }
+  
+  if (j < nlocal) {
+    for (m = 0; m < num_bond[j]; m++) {        
+    //  if(tag[i] == 25 or tag[j] == 25) printf("\t 2 ATOM LOOP %d-%d, %d/%d \n", tag[j], bond_atom[j][m], m, num_bond[j]);    
+      if (bond_atom[j][m] == tag[i]) {
+   //     if(tag[j] == 25 or tag[j] == 25) printf("\t DELETED\n");
+        bond_type[j][m] = 0;
+        n = num_bond[j];
+        bond_type[j][m] = bond_type[j][n-1];
+        bond_atom[j][m] = bond_atom[j][n-1];
+        num_bond[j]--;
+        break;
+      }      
+    }
+  }
+  
+ // if((tag[i] == 10 and tag[j] == 23)or (tag[i] == 23 and tag[j] == 10)) printf("Broken bond %d-%d (%d %d/%d)\n", tag[i], tag[j], i, j, nlocal);
+  
+}
+
+/* ----------------------------------------------------------------------
+   one method for every keyword bond bpm can output
+   the atom property is packed into array or vector
+------------------------------------------------------------------------- */
+
+void BondBPM::pack_id1(int n, int i, int j)
+{
+  tagint *tag = atom->tag;
+  output_data[n] = tag[i];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_id2(int n, int i, int j)
+{
+  tagint *tag = atom->tag;
+  output_data[n] = tag[j];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_time(int n, int i, int j)
+{
+  bigint time = update->ntimestep;
+  output_data[n] = time;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_x(int n, int i, int j)
+{
+  double **x = atom->x;
+  output_data[n] = (x[i][0] + x[j][0])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_y(int n, int i, int j)
+{
+  double **x = atom->x;
+  output_data[n] = (x[i][1] + x[j][1])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_z(int n, int i, int j)
+{
+  double **x = atom->x;
+  output_data[n] = (x[i][2] + x[j][2])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_x_ref(int n, int i, int j)
+{
+  double *x = atom->dvector[index_x_ref];
+  output_data[n] = (x[i] + x[j])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_y_ref(int n, int i, int j)
+{
+  double *y = atom->dvector[index_y_ref];
+  output_data[n] = (y[i] + y[j])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPM::pack_z_ref(int n, int i, int j)
+{
+  double *z = atom->dvector[index_z_ref];
+  output_data[n] = (z[i] + z[j])*0.5;
+}

src/BPM/bond_bpm.h

@@ -0,0 +1,90 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_BOND_BPM_H
+#define LMP_BOND_BPM_H
+
+#include "bond.h"
+
+namespace LAMMPS_NS {
+
+class BondBPM : public Bond {
+ public:
+  BondBPM(class LAMMPS *);
+  virtual ~BondBPM();
+  virtual void compute(int, int) = 0;
+  virtual void coeff(int, char **) = 0;
+  virtual void init_style();
+  void settings(int, char **);  
+  double equilibrium_distance(int);
+  void write_restart(FILE *){};
+  void read_restart(FILE *){};
+  void write_data(FILE *) {};
+  double single(int, double, int, int, double &) = 0;
+
+ protected:
+  double r0_max_estimate;
+  double max_stretch;
+  
+  char *id_fix_store_local, *id_fix_prop_atom;
+  class FixStoreLocal *fix_store_local;
+  class FixUpdateSpecialBonds *fix_update_special_bonds;
+  
+  void process_broken(int, int);  
+  typedef void (BondBPM::*FnPtrPack)(int,int,int);
+  FnPtrPack *pack_choice;              // ptrs to pack functions
+  double *output_data;
+
+  int prop_atom_flag, nvalues;
+  int index_x_ref, index_y_ref, index_z_ref;
+
+  void pack_id1(int,int,int);
+  void pack_id2(int,int,int);
+  void pack_time(int,int,int);
+  void pack_x(int,int,int);
+  void pack_y(int,int,int);
+  void pack_z(int,int,int);
+  void pack_x_ref(int,int,int);
+  void pack_y_ref(int,int,int);
+  void pack_z_ref(int,int,int);  
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Cannot find fix store/local
+
+Fix id cannot be found.
+
+E: Illegal bond_style command
+
+Self-explanatory.
+
+E: Bond style bpm/rotational must include at least one value to output
+
+Must include at least one bond property to store in fix store/local
+
+E: Bond style bpm/rotational cannot be used with 3,4-body interactions
+
+No angle, dihedral, or improper styles can be defined when using
+bond style bpm/rotational.
+
+E: Bond style bpm/rotational cannot be used with atom style template
+
+This bond style can change the bond topology which is not
+allowed with this atom style.
+
+*/

src/BPM/bond_bpm_rotational.cpp

@@ -0,0 +1,727 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "bond_bpm_rotational.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "error.h"
+#include "fix_bond_history.h"
+#include "force.h"
+#include "math_const.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "pair.h"
+
+#include <mpi.h>
+#include <cmath>
+#include <cstring>
+
+#define EPSILON 1e-10
+
+using namespace LAMMPS_NS;
+using namespace MathExtra;
+
+
+
+
+
+#include "update.h"
+/* ---------------------------------------------------------------------- */
+
+BondBPMRotational::BondBPMRotational(LAMMPS *lmp) : BondBPM(lmp)
+{
+  partial_flag = 1;
+  fix_bond_history = nullptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+BondBPMRotational::~BondBPMRotational()
+{
+  if(fix_bond_history) modify->delete_fix("BOND_HISTORY_BPM_ROTATIONAL");
+
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(Kr);
+    memory->destroy(Ks);
+    memory->destroy(Kt);
+    memory->destroy(Kb);
+    memory->destroy(Fcr);
+    memory->destroy(Fcs);
+    memory->destroy(Tct);
+    memory->destroy(Tcb);
+    memory->destroy(gnorm);
+    memory->destroy(gslide);
+    memory->destroy(groll);
+    memory->destroy(gtwist);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BondBPMRotational::acos_limit(double c)
+{
+  if (c > 1.0) c = 1.0;
+  if (c < -1.0) c = -1.0;
+  return acos(c);
+}
+
+/* ----------------------------------------------------------------------
+  Store data for a single bond - if bond added after LAMMPS init (e.g. pour)
+------------------------------------------------------------------------- */
+
+double BondBPMRotational::store_bond(int n,int i,int j)
+{
+  int m,k;
+  double delx, dely, delz, r, rinv;
+  double **x = atom->x;
+  tagint *tag = atom->tag;
+  double **bondstore = fix_bond_history->bondstore;
+
+  if (tag[i] < tag[j]) {
+    delx = x[i][0] - x[j][0];
+    dely = x[i][1] - x[j][1];
+    delz = x[i][2] - x[j][2];
+  } else {
+    delx = x[j][0] - x[i][0];
+    dely = x[j][1] - x[i][1];
+    delz = x[j][2] - x[i][2];
+  }
+
+  r = sqrt(delx*delx + dely*dely + delz*delz);
+  rinv = 1.0/r;
+
+  bondstore[n][0] = r;
+  bondstore[n][1] = delx*rinv;
+  bondstore[n][2] = dely*rinv;
+  bondstore[n][3] = delz*rinv;
+
+  if (i < atom->nlocal) {
+    for (m = 0; m < atom->num_bond[i]; m ++) {
+      if (atom->bond_atom[i][m] == tag[j]) {
+        fix_bond_history->update_atom_value(i, m, 0, r);
+        fix_bond_history->update_atom_value(i, m, 1, delx*rinv);
+        fix_bond_history->update_atom_value(i, m, 2, dely*rinv);
+        fix_bond_history->update_atom_value(i, m, 3, delz*rinv);
+      }
+    }
+  }
+
+  if (j < atom->nlocal) {
+    for (m = 0; m < atom->num_bond[j]; m ++) {
+      if (atom->bond_atom[j][m] == tag[i]) {
+        fix_bond_history->update_atom_value(j, m, 0, r);
+        fix_bond_history->update_atom_value(j, m, 1, delx*rinv);
+        fix_bond_history->update_atom_value(j, m, 2, dely*rinv);
+        fix_bond_history->update_atom_value(j, m, 3, delz*rinv);
+      }
+    }
+  }
+
+  return r;
+}
+
+/* ----------------------------------------------------------------------
+  Store data for all bonds called once
+------------------------------------------------------------------------- */
+
+void BondBPMRotational::store_data()
+{
+  int i, j, m, type;
+  double delx, dely, delz, r, rinv;
+  double **x = atom->x;
+  int **bond_type = atom->bond_type;
+  tagint *tag = atom->tag;
+
+  for (i = 0; i < atom->nlocal; i ++) {
+    for (m = 0; m < atom->num_bond[i]; m ++) {
+      type = bond_type[i][m];
+
+      //Skip if bond was turned off
+      if(type < 0)
+          continue;
+
+      // map to find index n for tag
+      j = atom->map(atom->bond_atom[i][m]);
+      if(j == -1) error->one(FLERR, "Atom missing in BPM bond");
+
+      // Save orientation as pointing towards small tag
+      if(tag[i] < tag[j]){
+        delx = x[i][0] - x[j][0];
+        dely = x[i][1] - x[j][1];
+        delz = x[i][2] - x[j][2];
+      } else {
+        delx = x[j][0] - x[i][0];
+        dely = x[j][1] - x[i][1];
+        delz = x[j][2] - x[i][2];
+      }
+      
+      // Get closest image in case bonded with ghost
+      domain->minimum_image(delx, dely, delz);
+      
+      r = sqrt(delx*delx + dely*dely + delz*delz);
+      rinv = 1.0/r;
+      fix_bond_history->update_atom_value(i, m, 0, r);
+      fix_bond_history->update_atom_value(i, m, 1, delx*rinv);
+      fix_bond_history->update_atom_value(i, m, 2, dely*rinv);
+      fix_bond_history->update_atom_value(i, m, 3, delz*rinv);
+    }
+  }
+
+  fix_bond_history->post_neighbor();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPMRotational::compute(int eflag, int vflag)
+{
+
+  if (! fix_bond_history->stored_flag) {
+    fix_bond_history->stored_flag = true;
+    store_data();
+  }
+
+  int i1,i2,itmp,m,n,type,itype,jtype;
+  double evdwl,fpair,rsq,ebond;
+  double q1[4], q2[4], r[3], r0[3];
+  double r0_mag, r_mag, r_mag_inv, Fr_mag, Fs_mag;
+  double Tt_mag, Tb_mag;
+  double force1on2[3], torque1on2[3], torque2on1[3];
+  double breaking, smooth, smooth_sq;
+  double rhat[3], wn1[3], wn2[3], wxn1[3], wxn2[3], vroll[3];
+  double w1dotr, w2dotr, v1dotr, v2dotr;
+  double vn1[3], vn2[3], vt1[3], vt2[3], tmp[3], s1[3], s2[3], tdamp[3]; 
+  double tor1, tor2, tor3, fs1, fs2, fs3;
+
+  double q2inv[4], rb[3], rb_x_r0[3], s[3], t[3], Fs[3];
+  double q21[4], qp21[4], Tbp[3], Ttp[3];
+  double Tsp[3], Fsp[3], Tt[3], Tb[3], Ts[3], F_rot[3], T_rot[3];
+  double mq[4], mqinv[4], Ttmp[3], Ftmp[3], qtmp[4];
+  double r0_dot_rb, gamma, c, psi, theta, sin_phi, cos_phi, temp, mag_in_plane, mag_out_plane;
+
+  ev_init(eflag,vflag);
+
+  if (vflag_global == 2)
+    force->pair->vflag_either = force->pair->vflag_global = 1;
+
+  double **cutsq = force->pair->cutsq;
+  double **x = atom->x;
+  double **v = atom->v;
+  double **omega = atom->omega;
+  double **f = atom->f;
+  double **torque = atom->torque;
+  double *radius = atom->radius;
+  double **quat = atom->quat;
+  tagint *tag = atom->tag;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+  int nlocal = atom->nlocal;
+  int newton_bond = force->newton_bond;
+
+  double **bondstore = fix_bond_history->bondstore;
+
+  for (n = 0; n < nbondlist; n++) {
+
+    // skip bond if already broken
+    if (bondlist[n][2] <= 0) continue;
+
+    i1 = bondlist[n][0];
+    i2 = bondlist[n][1];
+    type = bondlist[n][2];
+    r0_mag = bondstore[n][0];
+    
+    // Ensure pair is always ordered such that r0 points in
+    // a consistent direction and to ensure numerical operations
+    // are identical to minimize the possibility that a bond straddling
+    // an mpi grid (newton off) doesn't break on one proc but not the other
+    if (tag[i2] < tag[i1]) {
+      itmp = i1;
+      i1 = i2;
+      i2 = itmp;      
+    }
+
+    // If bond hasn't been set - should be initialized to zero
+    if (r0_mag < EPSILON || std::isnan(r0_mag))
+      r0_mag = store_bond(n,i1,i2);
+
+    r0[0] = bondstore[n][1];
+    r0[1] = bondstore[n][2];
+    r0[2] = bondstore[n][3];
+    MathExtra::scale3(r0_mag, r0);
+
+    q1[0] = quat[i1][0];
+    q1[1] = quat[i1][1];
+    q1[2] = quat[i1][2];
+    q1[3] = quat[i1][3];
+
+    q2[0] = quat[i2][0];
+    q2[1] = quat[i2][1];
+    q2[2] = quat[i2][2];
+    q2[3] = quat[i2][3];
+
+    // Note this is the reverse of Mora & Wang
+    MathExtra::sub3(x[i1], x[i2], r);
+
+    rsq = MathExtra::lensq3(r);
+    r_mag = sqrt(rsq);
+    r_mag_inv = 1.0/r_mag;
+    MathExtra::scale3(r_mag_inv, r, rhat);
+    
+    // ------------------------------------------------------//
+    //  Calculate forces using formulation in:
+    //    1) Y. Wang Acta Geotechnica 2009
+    //    2) P. Mora & Y. Wang Advances in Geomcomputing 2009
+    // ------------------------------------------------------//
+    
+    // Calculate normal forces, rb = bond vector in particle 1's frame
+    MathExtra::qconjugate(q2, q2inv);
+    MathExtra::quatrotvec(q2inv, r, rb);    
+    Fr_mag = Kr[type]*(r_mag - r0_mag);
+    
+    MathExtra::scale3(Fr_mag*r_mag_inv, rb, F_rot);
+    
+    // Calculate forces due to tangential displacements (no rotation)
+    r0_dot_rb = dot3(r0, rb);
+    c = r0_dot_rb*r_mag_inv/r0_mag;
+    gamma = acos_limit(c);
+    
+    MathExtra::cross3(rb, r0, rb_x_r0);
+    MathExtra::cross3(rb, rb_x_r0, s);
+    MathExtra::norm3(s);
+    
+    MathExtra::scale3(Ks[type]*r_mag*gamma, s, Fs);
+    
+    // Calculate torque due to tangential displacements
+    MathExtra::cross3(r0, rb, t);
+    MathExtra::norm3(t);
+    
+    MathExtra::scale3(0.5*r_mag*Ks[type]*r_mag*gamma, t, Ts);
+    
+    // Relative rotation force/torque
+    // Use representation of X'Y'Z' rotations from Wang, Mora 2009
+    temp = r_mag + rb[2];
+    if (temp < 0.0) temp = 0.0;
+    mq[0] = sqrt(2)*0.5*sqrt(temp*r_mag_inv);
+    
+    temp = sqrt(rb[0]*rb[0]+rb[1]*rb[1]);
+    if (temp != 0.0) {
+      mq[1] = -sqrt(2)*0.5/temp;
+      temp = r_mag - rb[2];
+      if (temp < 0.0) temp = 0.0;
+      mq[1] *= sqrt(temp*r_mag_inv);
+      mq[2] = -mq[1];
+      mq[1] *= rb[1];
+      mq[2] *= rb[0];
+    } else {
+      // If aligned along z axis, x,y terms zero (r_mag-rb[2] = 0)
+      mq[1] = 0.0;
+      mq[2] = 0.0;
+    }
+    mq[3] = 0.0;
+    
+    // qp21 = opposite of r^\circ_21 in Wang
+    // q21 = opposite of r_21 in Wang
+    MathExtra::quatquat(q2inv, q1, qp21);
+    MathExtra::qconjugate(mq, mqinv);
+    MathExtra::quatquat(mqinv,qp21,qtmp);
+    MathExtra::quatquat(qtmp,mq,q21);
+    
+    temp = sqrt(q21[0]*q21[0] + q21[3]*q21[3]);
+    if (temp != 0.0) {
+      c = q21[0]/temp;
+      psi = 2.0*acos_limit(c);
+    } else {
+      c = 0.0;
+      psi = 0.0;
+    }
+    
+    // Map negative rotations
+    if (q21[3] < 0.0) // sin = q21[3]/temp
+      psi = -psi;
+    
+    if (q21[3] == 0.0)
+      psi = 0.0;
+    
+    c = q21[0]*q21[0] - q21[1]*q21[1] - q21[2]*q21[2] + q21[3]*q21[3];
+    theta = acos_limit(c);
+    
+    // Separately calculte magnitude of quaternion in x-y and out of x-y planes
+    // to avoid dividing by zero
+    mag_out_plane = (q21[0]*q21[0] + q21[3]*q21[3]);
+    mag_in_plane = (q21[1]*q21[1] + q21[2]*q21[2]);
+    
+    if (mag_in_plane == 0.0) {
+      // No rotation => no bending/shear torque or extra shear force
+      // achieve by setting cos/sin = 0
+      cos_phi = 0.0;
+      sin_phi = 0.0;
+    } else if (mag_out_plane == 0.0) {
+      // Calculate angle in plane
+      cos_phi =  q21[2]/sqrt(mag_in_plane);
+      sin_phi = -q21[1]/sqrt(mag_in_plane);
+    } else {
+      // Default equations in Mora, Wang 2009
+      cos_phi = q21[1]*q21[3] + q21[0]*q21[2];
+      sin_phi = q21[2]*q21[3] - q21[0]*q21[1];
+    
+      cos_phi /= sqrt(mag_out_plane*mag_in_plane);
+      sin_phi /= sqrt(mag_out_plane*mag_in_plane);
+    }
+    
+    Tbp[0] = -Kb[type]*theta*sin_phi;
+    Tbp[1] = Kb[type]*theta*cos_phi;
+    Tbp[2] = 0.0;
+    
+    Ttp[0] = 0.0;
+    Ttp[1] = 0.0;
+    Ttp[2] = Kt[type]*psi;
+    
+    Fsp[0] = -0.5*Ks[type]*r_mag*theta*cos_phi;
+    Fsp[1] = -0.5*Ks[type]*r_mag*theta*sin_phi;
+    Fsp[2] = 0.0;
+    
+    Tsp[0] = 0.25*Ks[type]*r_mag*r_mag*theta*sin_phi;
+    Tsp[1] = -0.25*Ks[type]*r_mag*r_mag*theta*cos_phi;
+    Tsp[2] = 0.0;
+    
+    // Rotate forces/torques back to 1st particle's frame
+    
+    MathExtra::quatrotvec(mq, Fsp, Ftmp);
+    MathExtra::quatrotvec(mq, Tsp, Ttmp);
+    for (m = 0; m < 3; m++) {
+      Fs[m] += Ftmp[m];
+      Ts[m] += Ttmp[m];
+    }
+    
+    MathExtra::quatrotvec(mq, Tbp, Tb);
+    MathExtra::quatrotvec(mq, Ttp, Tt);
+    
+    // Sum forces and calculate magnitudes
+    F_rot[0] += Fs[0];
+    F_rot[1] += Fs[1];
+    F_rot[2] += Fs[2];
+    MathExtra::quatrotvec(q2, F_rot, force1on2);
+    
+    T_rot[0] = Ts[0] + Tt[0] + Tb[0];
+    T_rot[1] = Ts[1] + Tt[1] + Tb[1];
+    T_rot[2] = Ts[2] + Tt[2] + Tb[2];
+    MathExtra::quatrotvec(q2, T_rot, torque1on2);
+    
+    T_rot[0] = Ts[0] - Tt[0] - Tb[0];
+    T_rot[1] = Ts[1] - Tt[1] - Tb[1];
+    T_rot[2] = Ts[2] - Tt[2] - Tb[2];
+    MathExtra::quatrotvec(q2, T_rot, torque2on1);
+    
+    Fs_mag = MathExtra::len3(Fs);
+    Tt_mag = MathExtra::len3(Tt);
+    Tb_mag = MathExtra::len3(Tb);
+
+    // ------------------------------------------------------//
+    //  Check if bond breaks
+    // ------------------------------------------------------//
+    
+    if (r_mag < r0_mag)
+      breaking = Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type];
+    else
+      breaking = Fr_mag/Fcr[type] + Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type];
+
+    if (breaking >= 1.0) {        
+      bondlist[n][2] = 0;
+      process_broken(i1, i2);
+      continue;
+    }
+
+    smooth = 1.0 - breaking;
+    smooth_sq = smooth*smooth;
+
+    // Not actual energy, just a handy metric
+    if (eflag) ebond = -smooth_sq;
+
+    // ------------------------------------------------------//
+    //  Calculate damping using formulation in 
+    //    Y. Wang, F. Alonso-Marroquin, W. Guo 2015
+    // ------------------------------------------------------//
+    // Note: n points towards 1 vs pointing towards 2
+
+    // Damp normal velocity difference
+    v1dotr = MathExtra::dot3(v[i1],rhat);
+    v2dotr = MathExtra::dot3(v[i2],rhat);
+
+    MathExtra::scale3(v1dotr, rhat, vn1);
+    MathExtra::scale3(v2dotr, rhat, vn2);
+    
+    MathExtra::sub3(vn1, vn2, tmp);
+    MathExtra::scale3(gnorm[type], tmp);
+    MathExtra::add3(force1on2, tmp, force1on2);    
+    
+    // Damp tangential objective velocities    
+    MathExtra::sub3(v[i1], vn1, vt1);
+    MathExtra::sub3(v[i2], vn2, vt2);
+    
+    MathExtra::sub3(vt2, vt1, tmp);
+    MathExtra::scale3(-0.5, tmp);
+    
+    MathExtra::cross3(omega[i1], r, s1);
+    MathExtra::scale3(0.5, s1);
+    MathExtra::sub3(s1, tmp, s1); // Eq 12
+
+    MathExtra::cross3(omega[i2], r, s2);
+    MathExtra::scale3(-0.5,s2);
+    MathExtra::add3(s2, tmp, s2); // Eq 13
+    MathExtra::scale3(-0.5,s2);
+
+    MathExtra::sub3(s1, s2, tmp);
+    MathExtra::scale3(gslide[type], tmp);
+    MathExtra::add3(force1on2, tmp, force1on2);
+    
+    // Apply corresponding torque
+    MathExtra::cross3(r,tmp,tdamp);
+    MathExtra::scale3(-0.5, tdamp); // 0.5*r points from particle 2 to midpoint
+    MathExtra::add3(torque1on2, tdamp, torque1on2);
+    MathExtra::add3(torque2on1, tdamp, torque2on1);
+
+    // Damp rolling
+    MathExtra::cross3(omega[i1], rhat, wxn1);
+    MathExtra::cross3(omega[i2], rhat, wxn2);
+    MathExtra::sub3(wxn1, wxn2, vroll); // Eq. 31
+    MathExtra::cross3(r, vroll, tdamp); 
+
+    MathExtra::scale3(0.5*groll[type], tdamp);
+    MathExtra::add3(torque1on2, tdamp, torque1on2);
+    MathExtra::scale3(-1.0, tdamp);
+    MathExtra::add3(torque2on1, tdamp, torque2on1);
+
+    // Damp twist
+    w1dotr = MathExtra::dot3(omega[i1],rhat);
+    w2dotr = MathExtra::dot3(omega[i2],rhat);
+
+    MathExtra::scale3(w1dotr, rhat, wn1);
+    MathExtra::scale3(w2dotr, rhat, wn2);
+
+    MathExtra::sub3(wn1, wn2, tdamp); // Eq. 38
+    MathExtra::scale3(0.5*gtwist[type], tdamp);
+    MathExtra::add3(torque1on2, tdamp, torque1on2);
+    MathExtra::scale3(-1.0, tdamp);
+    MathExtra::add3(torque2on1, tdamp, torque2on1);
+
+    // ------------------------------------------------------//
+    //  Apply forces and torques to particles
+    // ------------------------------------------------------//
+
+    if (newton_bond || i1 < nlocal) {
+      f[i1][0] -= force1on2[0]*smooth_sq;
+      f[i1][1] -= force1on2[1]*smooth_sq;
+      f[i1][2] -= force1on2[2]*smooth_sq;
+
+      torque[i1][0] += torque2on1[0]*smooth_sq;
+      torque[i1][1] += torque2on1[1]*smooth_sq;
+      torque[i1][2] += torque2on1[2]*smooth_sq;
+    }
+
+    if (newton_bond || i2 < nlocal) {
+      f[i2][0] += force1on2[0]*smooth_sq;
+      f[i2][1] += force1on2[1]*smooth_sq;
+      f[i2][2] += force1on2[2]*smooth_sq;
+
+      torque[i2][0] += torque1on2[0]*smooth_sq;
+      torque[i2][1] += torque1on2[1]*smooth_sq;
+      torque[i2][2] += torque1on2[2]*smooth_sq;
+    }
+
+    if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,Fr_mag,r[0],r[1],r[2]);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPMRotational::allocate()
+{
+  allocated = 1;
+  int n = atom->nbondtypes;
+
+  memory->create(Kr,n+1,"bond:Kr");
+  memory->create(Ks,n+1,"bond:Ks");
+  memory->create(Kt,n+1,"bond:Kt");
+  memory->create(Kb,n+1,"bond:Kb");
+  memory->create(Fcr,n+1,"bond:Fcr");
+  memory->create(Fcs,n+1,"bond:Fcs");
+  memory->create(Tct,n+1,"bond:Tct");
+  memory->create(Tcb,n+1,"bond:Tcb");
+  memory->create(gnorm,n+1,"bond:gnorm");
+  memory->create(gslide,n+1,"bond:gslide");
+  memory->create(groll,n+1,"bond:groll");
+  memory->create(gtwist,n+1,"bond:gtwist");
+
+  memory->create(setflag,n+1,"bond:setflag");
+  for (int i = 1; i <= n; i++) setflag[i] = 0;
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more types
+------------------------------------------------------------------------- */
+
+void BondBPMRotational::coeff(int narg, char **arg)
+{
+  if (narg != 13) error->all(FLERR,"Incorrect args for bond coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi;
+  utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error);
+
+  double Kr_one = utils::numeric(FLERR,arg[1],false,lmp);
+  double Ks_one = utils::numeric(FLERR,arg[2],false,lmp);
+  double Kt_one = utils::numeric(FLERR,arg[3],false,lmp);
+  double Kb_one = utils::numeric(FLERR,arg[4],false,lmp);
+  double Fcr_one = utils::numeric(FLERR,arg[5],false,lmp);
+  double Fcs_one = utils::numeric(FLERR,arg[6],false,lmp);
+  double Tct_one = utils::numeric(FLERR,arg[7],false,lmp);
+  double Tcb_one = utils::numeric(FLERR,arg[8],false,lmp);
+  double gnorm_one = utils::numeric(FLERR,arg[9],false,lmp);
+  double gslide_one = utils::numeric(FLERR,arg[10],false,lmp);
+  double groll_one = utils::numeric(FLERR,arg[11],false,lmp);
+  double gtwist_one = utils::numeric(FLERR,arg[12],false,lmp);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    Kr[i] = Kr_one;
+    Ks[i] = Ks_one;
+    Kt[i] = Kt_one;
+    Kb[i] = Kb_one;
+    Fcr[i] = Fcr_one;
+    Fcs[i] = Fcs_one;
+    Tct[i] = Tct_one;
+    Tcb[i] = Tcb_one;
+    gnorm[i] = gnorm_one;
+    gslide[i] = gslide_one;
+    groll[i] = groll_one;
+    gtwist[i] = gtwist_one;
+    setflag[i] = 1;
+    count++;
+    
+    if (Fcr[i]/Kr[i] > max_stretch) max_stretch = Fcr[i]/Kr[i];
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   check if pair defined and special_bond settings are valid
+------------------------------------------------------------------------- */
+
+void BondBPMRotational::init_style()
+{
+  BondBPM::init_style();
+  
+  if (!atom->quat_flag || !atom->sphere_flag)
+    error->all(FLERR,"Bond bpm/rotational requires atom style sphere/bpm");
+  if (comm->ghost_velocity == 0)
+    error->all(FLERR,"Bond bpm/rotational requires ghost atoms store velocity");
+
+  if(domain->dimension == 2)
+    error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use");
+
+  if (!fix_bond_history)
+    fix_bond_history = (FixBondHistory *) modify->add_fix(
+         "BOND_HISTORY_BPM_ROTATIONAL all BOND_HISTORY 0 4");
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes out coeffs to restart file
+------------------------------------------------------------------------- */
+
+void BondBPMRotational::write_restart(FILE *fp)
+{
+  fwrite(&Kr[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Ks[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Kt[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Kb[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Fcr[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Fcs[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Tct[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Tcb[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&gnorm[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&gslide[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&groll[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&gtwist[1],sizeof(double),atom->nbondtypes,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads coeffs from restart file, bcasts them
+------------------------------------------------------------------------- */
+
+void BondBPMRotational::read_restart(FILE *fp)
+{
+  allocate();
+
+  if (comm->me == 0) {
+    utils::sfread(FLERR,&Kr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Ks[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Kt[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Kb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Fcr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Fcs[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Tct[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&Tcb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&gnorm[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&gslide[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&groll[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+    utils::sfread(FLERR,&gtwist[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
+  }
+  MPI_Bcast(&Kr[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Ks[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Kt[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Kb[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Fcr[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Fcs[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Tct[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Tcb[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&gnorm[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&gslide[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&groll[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&gtwist[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+
+  for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void BondBPMRotational::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nbondtypes; i++)
+     fprintf(fp,"%d %g %g %g %g %g %g %g %g %g %g %g %g\n",
+             i,Kr[i],Ks[i],Kt[i],Kb[i],Fcr[i], Fcs[i], Tct[i], 
+             Tcb[i], gnorm[i], gslide[i], groll[i], gtwist[i]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BondBPMRotational::single(int type, double rsq, int i, int j,
+                           double &fforce)
+{
+  // Not yet enabled
+  if (type <= 0) return 0.0;
+
+  //double r0;
+  //for (int n = 0; n < atom->num_bond[i]; n ++) {
+  //  if (atom->bond_atom[i][n] == atom->tag[j]) {
+  //    r0 = fix_bond_history->get_atom_value(i, n, 0);
+  //  }
+  //}
+}

src/BPM/bond_bpm_rotational.h

@@ -0,0 +1,77 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef BOND_CLASS
+// clang-format off
+BondStyle(bpm/rotational,BondBPMRotational)
+// clang-format on
+#else
+
+#ifndef LMP_BOND_BPM_ROTATIONAL_H
+#define LMP_BOND_BPM_ROTATIONAL_H
+
+#include "bond_bpm.h"
+
+namespace LAMMPS_NS {
+
+class BondBPMRotational : public BondBPM {
+ public:
+  BondBPMRotational(class LAMMPS *);
+  virtual ~BondBPMRotational();
+  virtual void compute(int, int);
+  void coeff(int, char **);
+  void init_style();
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_data(FILE *);
+  double single(int, double, int, int, double &);
+
+ protected:
+  double *Kr, *Ks, *Kt, *Kb, *gnorm, *gslide, *groll, *gtwist;
+  double *Fcr, *Fcs, *Tct, *Tcb;
+  double acos_limit(double);
+
+  void allocate();
+  void store_data();
+  double store_bond(int, int, int);
+  class FixBondHistory *fix_bond_history;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Atom missing in BPM bond
+
+Bonded atom cannot be found
+
+E: Incorrect args for bond coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Bond bpm/rotational requires atom style sphere/bpm
+
+Self-explanatory.
+
+E: Bond style bpm requires 1-3 and 1-4 special weights of 1.0
+
+Self-explanatory.
+
+W: Bond style bpm/rotational not intended for 2d use, may be inefficient
+
+This bond style will perform a lot of unnecessary calculations in 2d
+
+*/

src/BPM/compute_nbond_atom.cpp

@@ -0,0 +1,147 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_nbond_atom.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "update.h"
+
+#include <string.h>
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeNBondAtom::ComputeNBondAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  nbond(nullptr)
+{
+  if (narg < 3) error->all(FLERR,"Illegal compute nbond/atom command");
+
+  peratom_flag = 1;
+  size_peratom_cols = 0;
+  peatomflag = 1;
+  timeflag = 1;
+  comm_reverse = 1;
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeNBondAtom::~ComputeNBondAtom()
+{
+  memory->destroy(nbond);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeNBondAtom::compute_peratom()
+{
+
+  invoked_peratom = update->ntimestep;
+  if (update->eflag_atom != invoked_peratom)
+    error->all(FLERR,"Per-atom nbond was not tallied on needed timestep");
+
+  // grow local nbond array if necessary
+  // needs to be atom->nmax in length
+
+  if (atom->nmax > nmax) {
+    memory->destroy(nbond);
+    nmax = atom->nmax;
+    memory->create(nbond, nmax,"nbond/atom:nbond");
+    vector_atom = nbond;
+  }
+
+  // npair includes ghosts if either newton flag is set
+  //   b/c some bonds/dihedrals call pair::ev_tally with pairwise info
+  // nbond includes ghosts if newton_bond is set
+  // ntotal includes ghosts if either newton flag is set
+  // KSpace includes ghosts if tip4pflag is set
+
+  int nlocal = atom->nlocal;
+  tagint **bond_atom = atom->bond_atom;
+  int **bond_type = atom->bond_type;
+
+  int ntotal = nlocal;
+  if (force->newton) ntotal += atom->nghost;
+
+  // set local nbond array
+  int i, j, k;
+  int *num_bond = atom->num_bond;
+  int newton_bond = force->newton_bond;
+
+  for (i = 0; i < ntotal; i++) nbond[i] = 0;
+
+  for (i = 0; i < nlocal; i++) {
+    for (j = 0; j <num_bond[i]; j ++) {
+      if (bond_type[i][j] <= 0) continue;
+
+      k = atom->map(bond_atom[i][j]);
+      if (k < 0) continue;
+
+      nbond[i] += 1;
+      if (newton_bond) nbond[k] += 1;
+    }
+  }
+
+  // communicate ghost nbond between neighbor procs
+  if (force->newton)
+    comm->reverse_comm_compute(this);
+
+  // zero nbond of atoms not in group
+  // only do this after comm since ghost contributions must be included
+  int *mask = atom->mask;
+
+  for (i = 0; i < nlocal; i++)
+    if (!(mask[i] & groupbit)) nbond[i] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeNBondAtom::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) buf[m++] = nbond[i];
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeNBondAtom::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    nbond[j] += buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeNBondAtom::memory_usage()
+{
+  double bytes = nmax * sizeof(double);
+  return bytes;
+}

src/BPM/compute_nbond_atom.h

@@ -0,0 +1,61 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(nbond/atom,ComputeNBondAtom)
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_NBOND_ATOM_H
+#define LMP_COMPUTE_NBOND_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeNBondAtom : public Compute {
+ public:
+  ComputeNBondAtom(class LAMMPS *, int, char **);
+  ~ComputeNBondAtom();
+  void init() {}
+  void compute_peratom();
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  double memory_usage();
+
+ private:
+  int nmax;
+  double *nbond;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Per-atom energy was not tallied on needed timestep
+
+You are using a thermo keyword that requires potentials to
+have tallied energy, but they didn't on this timestep.  See the
+variable doc page for ideas on how to make this work.
+
+*/

src/BPM/fix_bond_history.cpp

@@ -0,0 +1,294 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "fix_bond_history.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "memory.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "string.h"
+
+#include <cstring>
+#include <mpi.h>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define LB_FACTOR 1.5
+#define DELTA 10000
+
+/* ---------------------------------------------------------------------- */
+
+FixBondHistory::FixBondHistory(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 5) error->all(FLERR,"Illegal fix bond/history command");
+  update_flag = utils::inumeric(FLERR,arg[3],false,lmp);
+  ndata = utils::inumeric(FLERR,arg[4],false,lmp);
+  nbond = atom->bond_per_atom;
+
+  if (nbond == 0)
+    error->all(FLERR, "Cannot store bond variables without any bonds");
+
+  stored_flag = false;
+  restart_global = 1;
+  create_attribute = 1;
+
+  bondstore = nullptr;
+  maxbond = 0;
+  allocate();
+
+  new_fix_id = nullptr;
+  array_id = nullptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixBondHistory::~FixBondHistory()
+{
+  if (new_fix_id && modify->nfix) modify->delete_fix(new_fix_id);
+  delete [] new_fix_id;
+  delete [] array_id;
+
+  memory->destroy(bondstore);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixBondHistory::setmask()
+{
+  int mask = 0;
+  mask |= PRE_EXCHANGE;
+  mask |= POST_NEIGHBOR;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondHistory::post_constructor()
+{
+  // Store saved bond quantities for each atom using fix property atom
+  char **newarg = new char*[5];
+  int nvalue = 0;
+  int tmp1, tmp2;
+
+  int nn = strlen(id) + strlen("_FIX_PROP_ATOM") + 1;
+  new_fix_id = new char[nn];
+  strcpy(new_fix_id, "FIX_PROP_ATOM_");
+  strcat(new_fix_id, id);
+
+  nn = strlen(id) + 4 ;
+  array_id = new char[nn];
+  strcpy(array_id, "d2_");
+  strcat(array_id, id);
+
+  char ncols[16];
+  sprintf(ncols,"%d",nbond*ndata);
+  newarg[0] = new_fix_id;
+  newarg[1] = group->names[igroup];
+  newarg[2] = (char *) "property/atom";
+  newarg[3] = array_id;
+  newarg[4] = ncols;
+
+  modify->add_fix(5,newarg);
+  index = atom->find_custom(&array_id[3],tmp1,tmp2);
+
+  delete [] newarg;
+} 
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondHistory::update_atom_value(int i, int m, int idata, double value)
+{
+  if (idata >= ndata || m > nbond) error->all(FLERR, "Index exceeded in fix bond history");
+  atom->darray[index][i][m*ndata+idata] = value;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixBondHistory::get_atom_value(int i, int m, int idata)
+{
+   if (idata >= ndata || m > nbond) error->all(FLERR, "Index exceeded in fix bond history");
+   return atom->darray[index][i][m*ndata+idata];
+}
+
+/* ----------------------------------------------------------------------
+   If stored values are updated, need to copy to atom arrays before exchanging
+   Always called before neighborlist rebuilt
+   Also call prior to irregular communication in other fixes (e.g. deform)
+------------------------------------------------------------------------- */
+
+void FixBondHistory::pre_exchange()
+{
+  if (!update_flag) return;
+
+  int i1, i2, n, m, idata;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+  double **stored = atom->darray[index];
+
+  int nlocal = atom->nlocal;
+  tagint **bond_atom = atom->bond_atom;
+  int *num_bond = atom->num_bond;
+  tagint *tag = atom->tag;
+
+  for (n = 0; n < nbondlist; n++) {
+    i1 = bondlist[n][0];
+    i2 = bondlist[n][1];
+
+    // skip bond if already broken
+    if (bondlist[n][2] <= 0) {
+      continue;
+    }
+
+    if (i1 < nlocal) {
+      for (m = 0; m < num_bond[i1]; m++) {
+        if (bond_atom[i1][m] == tag[i2]) {
+          for (idata = 0; idata < ndata; idata++) {
+            stored[i1][m*ndata+idata] = bondstore[n][idata];
+          }
+        }
+      }
+    }
+
+    if (i2 < nlocal) {
+      for (m = 0; m < num_bond[i2]; m++) {
+        if (bond_atom[i2][m] == tag[i1]) {
+          for (idata = 0; idata < ndata; idata++) {
+            stored[i1][m*ndata+idata] = bondstore[n][idata];
+          }
+        }
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondHistory::allocate()
+{
+  //Ideally would just ask ntopo for maxbond, but protected
+  if (comm->nprocs == 1) maxbond = atom->nbonds;
+  else maxbond = static_cast<int> (LB_FACTOR * atom->nbonds / comm->nprocs);
+  memory->create(bondstore,maxbond,ndata,"fix_bond_store:bondstore");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondHistory::setup_post_neighbor()
+{
+  post_neighbor();
+}
+
+/* ----------------------------------------------------------------------
+   called after neighbor list is build
+   build array of stored bond quantities from fix property atom
+------------------------------------------------------------------------- */
+
+void FixBondHistory::post_neighbor()
+{
+  //Grow array if number of bonds has increased
+  while (neighbor->nbondlist >= maxbond) {
+    maxbond += DELTA;
+    memory->grow(bondstore,maxbond,ndata,"fix_bond_store:bondstore");
+  }
+
+  int i1, i2, n, m, idata;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+  double **stored = atom->darray[index];
+
+  int nlocal = atom->nlocal;
+  tagint **bond_atom = atom->bond_atom;
+  int *num_bond = atom->num_bond;
+  tagint *tag = atom->tag;
+
+  for (n = 0; n < nbondlist; n++) {
+    i1 = bondlist[n][0];
+    i2 = bondlist[n][1];
+
+    // skip bond if already broken
+    if (bondlist[n][2] <= 0) {
+      continue;
+    }
+
+    if (i1 < nlocal) {
+      for (m = 0; m < num_bond[i1]; m++) {
+        if (bond_atom[i1][m] == tag[i2]) {
+          for (idata = 0; idata < ndata; idata++) {
+            bondstore[n][idata] = stored[i1][m*ndata+idata];
+          }
+        }
+      }
+    }
+
+    if (i2 < nlocal){
+      for (m = 0; m < num_bond[i2]; m++) {
+        if (bond_atom[i2][m] == tag[i1]) {
+          for (idata = 0; idata < ndata; idata++) {
+            bondstore[n][idata] = stored[i2][m*ndata+idata];
+          }
+        }
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixBondHistory::memory_usage()
+{
+  return maxbond * ndata * sizeof(double);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondHistory::write_restart(FILE *fp)
+{
+  int n = 0;
+  double list[1];
+  list[n++] = stored_flag;
+
+  if (comm->me == 0) {
+    int size = n * sizeof(double);
+    fwrite(&size,sizeof(int),1,fp);
+    fwrite(list,sizeof(double),n,fp);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondHistory::restart(char *buf)
+{
+  int n = 0;
+  double *list = (double *) buf;
+  stored_flag = static_cast<int> (list[n++]);
+}
+
+/* ----------------------------------------------------------------------
+   initialize bond values to zero, called when atom is created
+------------------------------------------------------------------------- */
+
+void FixBondHistory::set_arrays(int i)
+{
+  double **stored = atom->darray[index];
+  for (int m = 0; m < nbond; m++)
+    for (int idata = 0; idata < ndata; idata++)
+      stored[i][m*ndata+idata] = 0.0;
+}

src/BPM/fix_bond_history.h

@@ -0,0 +1,80 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(BOND_HISTORY,FixBondHistory)
+// clang-format on
+#else
+
+#ifndef LMP_FIX_BOND_HISTORY_H
+#define LMP_FIX_BOND_HISTORY_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixBondHistory : public Fix {
+ public:
+
+  FixBondHistory(class LAMMPS *, int, char **);
+  ~FixBondHistory();
+  int setmask();
+  void post_constructor();
+  void setup_post_neighbor();
+  void post_neighbor();
+  void pre_exchange();
+  double memory_usage();
+  void write_restart(FILE *fp);
+  void restart(char *buf);
+  void set_arrays(int);
+
+  void update_atom_value(int, int, int, double);
+  double get_atom_value(int, int, int);
+  double **bondstore;
+  int stored_flag;
+
+ protected:
+
+  void allocate();
+
+  int update_flag;
+  int nbond, maxbond, ndata;
+  int index;
+  char *new_fix_id;
+  char *array_id;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+UNDOCUMENTED
+
+E: Index exceeded in fix bond history
+
+Bond requested non-existant data
+
+E: Cannot store bond variables without any bonds
+
+Atoms must have a nonzero number of bonds to store data
+
+
+
+*/
+

src/BPM/fix_nve_sphere_bpm.cpp

@@ -0,0 +1,161 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "fix_nve_sphere_bpm.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "math_extra.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathExtra;
+
+
+/* ---------------------------------------------------------------------- */
+
+FixNVESphereBPM::FixNVESphereBPM(LAMMPS *lmp, int narg, char **arg) :
+  FixNVE(lmp, narg, arg)
+{
+  if (narg < 3) error->all(FLERR,"Illegal fix nve/sphere/bpm command");
+
+  time_integrate = 1;
+
+  // process extra keywords
+  // inertia = moment of inertia prefactor for sphere or disc
+
+  inertia = 0.4;
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"disc")==0) {
+      inertia = 0.5;
+      if (domain->dimension != 2)
+        error->all(FLERR,"Fix nve/sphere/bpm disc requires 2d simulation");
+      iarg++;
+    }
+    else error->all(FLERR,"Illegal fix nve/sphere/bpm command");
+  }
+
+  inv_inertia = 1.0/inertia;
+
+  // error checks
+
+  if (!atom->quat_flag || !atom->sphere_flag)
+    error->all(FLERR,"Fix nve/sphere/bpm requires atom style sphere/bpm");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVESphereBPM::init()
+{
+  FixNVE::init();
+
+  // check that all particles are finite-size spheres
+  // no point particles allowed
+
+  double *radius = atom->radius;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit)
+      if (radius[i] == 0.0)
+        error->one(FLERR,"Fix nve/sphere/bpm requires extended particles");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVESphereBPM::initial_integrate(int /*vflag*/)
+{
+  double dtq,dtfm,dtirotate,particle_inertia;
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **omega = atom->omega;
+  double **torque = atom->torque;
+  double **quat = atom->quat;
+  double *radius = atom->radius;
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // set timestep here since dt may have changed or come via rRESPA
+  dtq = 0.5 * dtv;
+
+  // update v,x,omega,quat for all particles
+  // d_omega/dt = torque / inertia
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / rmass[i];
+      v[i][0] += dtfm * f[i][0];
+      v[i][1] += dtfm * f[i][1];
+      v[i][2] += dtfm * f[i][2];
+      x[i][0] += dtv * v[i][0];
+      x[i][1] += dtv * v[i][1];
+      x[i][2] += dtv * v[i][2];
+
+      particle_inertia = inertia*(radius[i]*radius[i]*rmass[i]);
+      dtirotate = dtf / particle_inertia;
+      omega[i][0] += dtirotate * torque[i][0];
+      omega[i][1] += dtirotate * torque[i][1];
+      omega[i][2] += dtirotate * torque[i][2];
+
+      MathExtra::richardson_sphere(quat[i],omega[i],dtq);
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVESphereBPM::final_integrate()
+{
+  double dtfm,dtirotate,particle_inertia;
+
+  double **v = atom->v;
+  double **f = atom->f;
+  double **omega = atom->omega;
+  double **torque = atom->torque;
+  double *rmass = atom->rmass;
+  double *radius = atom->radius;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // update v,omega for all particles
+  // d_omega/dt = torque / inertia
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      dtfm = dtf / rmass[i];
+      v[i][0] += dtfm * f[i][0];
+      v[i][1] += dtfm * f[i][1];
+      v[i][2] += dtfm * f[i][2];
+
+      particle_inertia = inertia*(radius[i]*radius[i]*rmass[i]);
+      dtirotate = dtf / particle_inertia;
+      omega[i][0] += dtirotate * torque[i][0];
+      omega[i][1] += dtirotate * torque[i][1];
+      omega[i][2] += dtirotate * torque[i][2];
+    }
+}

src/BPM/fix_nve_sphere_bpm.h

@@ -0,0 +1,71 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(nve/sphere/bpm,FixNVESphereBPM)
+// clang-format on
+#else
+
+#ifndef LMP_FIX_NVE_SPHERE_BPM_H
+#define LMP_FIX_NVE_SPHERE_BPM_H
+
+#include "fix_nve.h"
+
+namespace LAMMPS_NS {
+
+class FixNVESphereBPM : public FixNVE {
+ public:
+  FixNVESphereBPM(class LAMMPS *, int, char **);
+  virtual ~FixNVESphereBPM() {}
+  void init();
+  virtual void initial_integrate(int);
+  virtual void final_integrate();
+
+ protected:
+  double inertia, inv_inertia;
+  int extra;
+  int dlm;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix nve/sphere/bpm disc requires 2d simulation
+
+UNDOCUMENTED
+
+E: Fix nve/sphere/bpm requires atom style sphere/bpm
+
+Self-explanatory.
+
+E: Fix nve/sphere/bpm update dipole requires atom attribute mu
+
+An atom style with this attribute is needed.
+
+E: Fix nve/sphere/bpm requires extended particles
+
+This fix can only be used for particles of a finite size.
+
+
+*/

src/BPM/fix_update_special_bonds.cpp

@@ -0,0 +1,235 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://lammps.sandia.gov/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "fix_update_special_bonds.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "pair.h"
+
+#include <cstring>
+#include <set>
+#include <utility>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixUpdateSpecialBonds::FixUpdateSpecialBonds(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal fix update/special/bonds command");
+  comm_forward = 1+atom->maxspecial;  
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixUpdateSpecialBonds::~FixUpdateSpecialBonds()
+{
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixUpdateSpecialBonds::setmask()
+{
+  int mask = 0;
+  mask |= PRE_EXCHANGE;
+  mask |= PRE_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixUpdateSpecialBonds::setup(int /*vflag*/)
+{
+  // Require atoms know about all of their bonds and if they break
+  if (force->newton_bond)
+    error->all(FLERR,"Fix update/special/bonds requires Newton bond off");
+
+  if (!atom->avec->bonds_allow)
+    error->all(FLERR,"Fix update/special/bonds requires atom bonds");
+
+  // special lj must be 0 1 1 to censor pair forces between bonded particles
+  // special coulomb must be 1 1 1 to ensure all pairs are included in the 
+  //   neighbor list and 1-3 and 1-4 special bond lists are skipped  
+  if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
+      force->special_lj[3] != 1.0)
+    error->all(FLERR,"Fix update/special/bonds requires special LJ weights = 0,1,1");
+  if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 ||
+      force->special_coul[3] != 1.0)
+    error->all(FLERR,"Fix update/special/bonds requires special Coulomb weights = 1,1,1");
+    
+  broken_pairs.clear();
+}
+
+/* ----------------------------------------------------------------------
+  Update special bond list and atom bond arrays, empty broken bond list
+------------------------------------------------------------------------- */
+
+void FixUpdateSpecialBonds::pre_exchange()
+{
+  int i, j, key, m, n1, n3;
+  tagint min_tag, max_tag;
+  int nlocal = atom->nlocal;
+  
+  tagint *tag = atom->tag;  
+  tagint *slist;
+  int **nspecial = atom->nspecial;
+  tagint **special = atom->special;
+  
+  for (auto const &key : broken_pairs) {
+    min_tag = key.first;
+    max_tag = key.second;
+    
+    i = atom->map(min_tag);
+    j = atom->map(max_tag);
+    
+    // remove i from special bond list for atom j and vice versa
+    if (i < nlocal) {
+      slist = special[i];
+      n1 = nspecial[i][0];
+      for (m = 0; m < n1; m++)
+        if (slist[m] == max_tag) break;
+      n3 = nspecial[i][2];
+      for (; m < n3-1; m++) slist[m] = slist[m+1];
+      nspecial[i][0]--;
+      nspecial[i][1]--;
+      nspecial[i][2]--;
+    }
+    
+    if (j < nlocal) {
+      slist = special[j];
+      n1 = nspecial[j][0];
+      for (m = 0; m < n1; m++)
+        if (slist[m] == min_tag) break;
+      n3 = nspecial[j][2];
+      for (; m < n3-1; m++) slist[m] = slist[m+1];
+      nspecial[j][0]--;
+      nspecial[j][1]--;
+      nspecial[j][2]--;
+    }
+  }
+  
+  // Forward updated special bond list
+  comm->forward_comm_fix(this);
+  
+  broken_pairs.clear();
+}
+
+/* ----------------------------------------------------------------------
+  Loop neighbor list and update special bond lists for recently broken bonds
+------------------------------------------------------------------------- */
+
+void FixUpdateSpecialBonds::pre_force(int /*vflag*/)
+{
+  int i,j,n,m,ii,jj,inum,jnum;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  tagint min_tag, max_tag;
+  std::pair <tagint, tagint> key;
+  
+  int **bond_type = atom->bond_type;
+  int *num_bond = atom->num_bond;
+  tagint **bond_atom = atom->bond_atom;
+  int nlocal = atom->nlocal;
+  
+  tagint *tag = atom->tag;
+  NeighList *list = force->pair->list;  
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      
+      min_tag = tag[i];
+      max_tag = tag[j];
+      if (max_tag < min_tag) {
+        min_tag = tag[j];
+        max_tag = tag[i];
+      }
+      key = std::make_pair(min_tag, max_tag);
+
+      if (broken_pairs.find(key) != broken_pairs.end()) 
+        jlist[jj] = j; // Clear special bond bits
+    } 
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixUpdateSpecialBonds::pack_forward_comm(int n, int *list, double *buf,
+                                    int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i,j,k,m,ns;
+  int **nspecial = atom->nspecial;
+  tagint **special = atom->special;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    ns = nspecial[j][0];
+    buf[m++] = ubuf(ns).d;
+    for (k = 0; k < ns; k++)
+      buf[m++] = ubuf(special[j][k]).d;
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixUpdateSpecialBonds::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i,j,m,ns,last;
+  int **nspecial = atom->nspecial;
+  tagint **special = atom->special;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    ns = (int) ubuf(buf[m++]).i;
+    nspecial[i][0] = ns;
+    for (j = 0; j < ns; j++)
+      special[i][j] = (tagint) ubuf(buf[m++]).i;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixUpdateSpecialBonds::add_broken_bond(int i, int j)
+{
+  tagint *tag = atom->tag;
+  
+  tagint min_tag = tag[i];
+  tagint max_tag = tag[j];
+  if (max_tag < min_tag) {
+    min_tag = tag[j];
+    max_tag = tag[i];
+  }
+  std::pair <tagint, tagint> key = std::make_pair(min_tag, max_tag);
+
+  broken_pairs.insert(key);
+}

src/BPM/fix_update_special_bonds.h

@@ -0,0 +1,76 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(update/special/bonds,FixUpdateSpecialBonds)
+// clang-format on
+#else
+
+#ifndef LMP_FIX_UPDATE_SPECIAL_BONDS_H
+#define LMP_FIX_UPDATE_SPECIAL_BONDS_H
+
+#include "fix.h"
+
+#include <set>
+#include <utility>
+
+namespace LAMMPS_NS {
+
+class FixUpdateSpecialBonds : public Fix {
+ public:
+  FixUpdateSpecialBonds(class LAMMPS *, int, char **);
+  ~FixUpdateSpecialBonds();
+  int setmask();
+  void setup(int);
+  void pre_exchange();
+  void pre_force(int);
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);  
+  void add_broken_bond(int,int);
+
+ protected:
+  std::set <std::pair<tagint, tagint>> broken_pairs;
+  inline int sbmask(int j) const {
+    return j >> SBBITS & 3;
+  }
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal fix censor/bonded/pairs command
+
+Self-explanatory.
+
+E: Fix censor/bonded/pairs requires Newton bond off
+
+Self-explanatory.
+
+E: Fix censor/bonded/pairs requires atom bonds
+
+Self-explanatory.
+
+E: Fix censor/bonded/pairs must be used without special bonds
+
+Self-explanatory. Look at the atom modify special command.
+
+E: Fix censor/bonded/pairs requires special_bonds = 0,0,0
+
+Self-explanatory.
+
+*/

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -43,7 +43,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
 {
   int i,j,ii,jj,inum,jnum;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
-  double radi,radj,radsum,rsq,r,rinv,rsqinv;
+  double radi,radj,radsum,rsq,r,rinv,rsqinv,factor_lj;
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
   double wr1,wr2,wr3;
   double vtr1,vtr2,vtr3,vrel;
@@ -112,8 +112,11 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
       j &= NEIGHMASK;
 
+      if (factor_lj == 0) continue;
+
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
@@ -244,6 +247,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
         fx = delx*ccel + fs1;
         fy = dely*ccel + fs2;
         fz = delz*ccel + fs3;
+        fx *= factor_lj;
+        fy *= factor_lj;
+        fz *= factor_lj;
         f[i][0] += fx;
         f[i][1] += fy;
         f[i][2] += fz;
@@ -251,6 +257,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
         tor1 = rinv * (dely*fs3 - delz*fs2);
         tor2 = rinv * (delz*fs1 - delx*fs3);
         tor3 = rinv * (delx*fs2 - dely*fs1);
+        tor1 *= factor_lj;
+        tor2 *= factor_lj;
+        tor3 *= factor_lj;
         torque[i][0] -= radi*tor1;
         torque[i][1] -= radi*tor2;
         torque[i][2] -= radi*tor3;

src/GRANULAR/pair_gran_hooke.cpp

@@ -42,7 +42,7 @@ void PairGranHooke::compute(int eflag, int vflag)
 {
   int i,j,ii,jj,inum,jnum;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
-  double radi,radj,radsum,rsq,r,rinv,rsqinv;
+  double radi,radj,radsum,rsq,r,rinv,rsqinv,factor_lj;
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
   double wr1,wr2,wr3;
   double vtr1,vtr2,vtr3,vrel;
@@ -101,8 +101,11 @@ void PairGranHooke::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
       j &= NEIGHMASK;
 
+      if (factor_lj == 0) continue;
+
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
@@ -187,6 +190,9 @@ void PairGranHooke::compute(int eflag, int vflag)
         fx = delx*ccel + fs1;
         fy = dely*ccel + fs2;
         fz = delz*ccel + fs3;
+        fx *= factor_lj;
+        fy *= factor_lj;
+        fz *= factor_lj;
         f[i][0] += fx;
         f[i][1] += fy;
         f[i][2] += fz;
@@ -194,6 +200,9 @@ void PairGranHooke::compute(int eflag, int vflag)
         tor1 = rinv * (dely*fs3 - delz*fs2);
         tor2 = rinv * (delz*fs1 - delx*fs3);
         tor3 = rinv * (delx*fs2 - dely*fs1);
+        tor1 *= factor_lj;
+        tor2 *= factor_lj;
+        tor3 *= factor_lj;
         torque[i][0] -= radi*tor1;
         torque[i][1] -= radi*tor2;
         torque[i][2] -= radi*tor3;

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -101,7 +101,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 {
   int i,j,ii,jj,inum,jnum;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
-  double radi,radj,radsum,rsq,r,rinv,rsqinv;
+  double radi,radj,radsum,rsq,r,rinv,rsqinv,factor_lj;
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
   double wr1,wr2,wr3;
   double vtr1,vtr2,vtr3,vrel;
@@ -170,8 +170,11 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
       j &= NEIGHMASK;
 
+      if (factor_lj == 0) continue;
+
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
@@ -301,6 +304,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
         fx = delx*ccel + fs1;
         fy = dely*ccel + fs2;
         fz = delz*ccel + fs3;
+        fx *= factor_lj;
+        fy *= factor_lj;
+        fz *= factor_lj;
         f[i][0] += fx;
         f[i][1] += fy;
         f[i][2] += fz;
@@ -308,6 +314,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
         tor1 = rinv * (dely*fs3 - delz*fs2);
         tor2 = rinv * (delz*fs1 - delx*fs3);
         tor3 = rinv * (delx*fs2 - dely*fs1);
+        tor1 *= factor_lj;
+        tor2 *= factor_lj;
+        tor3 *= factor_lj;
         torque[i][0] -= radi*tor1;
         torque[i][1] -= radi*tor2;
         torque[i][2] -= radi*tor3;

src/GRANULAR/pair_granular.cpp

@@ -150,7 +150,7 @@ void PairGranular::compute(int eflag, int vflag)
 {
   int i,j,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz;
-  double radi,radj,radsum,rsq,r,rinv;
+  double radi,radj,radsum,rsq,r,rinv,factor_lj;
   double Reff, delta, dR, dR2, dist_to_contact;
 
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -249,8 +249,11 @@ void PairGranular::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
       j &= NEIGHMASK;
 
+      if (factor_lj == 0) continue;
+
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
@@ -651,6 +654,9 @@ void PairGranular::compute(int eflag, int vflag)
         fx = nx*Fntot + fs1;
         fy = ny*Fntot + fs2;
         fz = nz*Fntot + fs3;
+        fx *= factor_lj;
+        fy *= factor_lj;
+        fz *= factor_lj;
 
         f[i][0] += fx;
         f[i][1] += fy;
@@ -659,6 +665,9 @@ void PairGranular::compute(int eflag, int vflag)
         tor1 = ny*fs3 - nz*fs2;
         tor2 = nz*fs1 - nx*fs3;
         tor3 = nx*fs2 - ny*fs1;
+        tor1 *= factor_lj;
+        tor2 *= factor_lj;
+        tor3 *= factor_lj;
 
         dist_to_contact = radi-0.5*delta;
         torque[i][0] -= dist_to_contact*tor1;
@@ -666,9 +675,9 @@ void PairGranular::compute(int eflag, int vflag)
         torque[i][2] -= dist_to_contact*tor3;
 
         if (twist_model[itype][jtype] != TWIST_NONE) {
-          tortwist1 = magtortwist * nx;
-          tortwist2 = magtortwist * ny;
-          tortwist3 = magtortwist * nz;
+          tortwist1 = magtortwist * nx * factor_lj;
+          tortwist2 = magtortwist * ny * factor_lj;
+          tortwist3 = magtortwist * nz * factor_lj;
 
           torque[i][0] += tortwist1;
           torque[i][1] += tortwist2;
@@ -676,9 +685,9 @@ void PairGranular::compute(int eflag, int vflag)
         }
 
         if (roll_model[itype][jtype] != ROLL_NONE) {
-          torroll1 = Reff*(ny*fr3 - nz*fr2); // n cross fr
-          torroll2 = Reff*(nz*fr1 - nx*fr3);
-          torroll3 = Reff*(nx*fr2 - ny*fr1);
+          torroll1 = Reff*(ny*fr3 - nz*fr2) * factor_lj; // n cross fr
+          torroll2 = Reff*(nz*fr1 - nx*fr3) * factor_lj;
+          torroll3 = Reff*(nx*fr2 - ny*fr1) * factor_lj;
 
           torque[i][0] += torroll1;
           torque[i][1] += torroll2;

src/MISC/fix_pair_tracker.cpp

@@ -1,369 +0,0 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   https://www.lammps.org/ Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#include "fix_pair_tracker.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "error.h"
-#include "group.h"
-#include "memory.h"
-#include "modify.h"
-#include "tokenizer.h"
-#include "update.h"
-
-#include <string.h>
-
-using namespace LAMMPS_NS;
-using namespace FixConst;
-
-#define DELTA 1000
-
-/* ---------------------------------------------------------------------- */
-
-FixPairTracker::FixPairTracker(LAMMPS *lmp, int narg, char **arg) :
-    Fix(lmp, narg, arg), nvalues(0), vector(NULL), array(NULL), pack_choice(NULL)
-{
-  if (narg < 3) error->all(FLERR, "Illegal fix pair/tracker command");
-  local_flag = 1;
-
-  nevery = utils::inumeric(FLERR, arg[3], false, lmp);
-  if (nevery <= 0) error->all(FLERR, "Illegal fix pair/tracker command");
-  local_freq = nevery;
-
-  // If optional arguments included, this will be oversized
-  nvalues = narg - 4;
-  pack_choice = new FnPtrPack[nvalues];
-
-  tmin = -1;
-  type_filter = nullptr;
-  int iarg = 4;
-  nvalues = 0;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg], "id1") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_id1;
-    } else if (strcmp(arg[iarg], "id2") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_id2;
-
-    } else if (strcmp(arg[iarg], "time/created") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_time_created;
-    } else if (strcmp(arg[iarg], "time/broken") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_time_broken;
-    } else if (strcmp(arg[iarg], "time/total") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_time_total;
-
-    } else if (strcmp(arg[iarg], "x") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_x;
-    } else if (strcmp(arg[iarg], "y") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_y;
-    } else if (strcmp(arg[iarg], "z") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_z;
-
-    } else if (strcmp(arg[iarg], "r/min") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_rmin;
-    } else if (strcmp(arg[iarg], "r/ave") == 0) {
-      pack_choice[nvalues++] = &FixPairTracker::pack_rave;
-
-    } else if (strcmp(arg[iarg], "time/min") == 0) {
-      if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
-      tmin = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      iarg++;
-
-    } else if (strcmp(arg[iarg], "type/include") == 0) {
-      if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
-      int ntypes = atom->ntypes;
-
-      int i, j, itype, jtype, in, jn, infield, jnfield;
-      int inlo, inhi, jnlo, jnhi;
-      char *istr, *jstr;
-      if (!type_filter) {
-        memory->create(type_filter, ntypes + 1, ntypes + 1, "fix/pair/tracker:type_filter");
-
-        for (i = 0; i <= ntypes; i++) {
-          for (j = 0; j <= ntypes; j++) { type_filter[i][j] = 0; }
-        }
-      }
-
-      in = strlen(arg[iarg + 1]) + 1;
-      istr = new char[in];
-      strcpy(istr, arg[iarg + 1]);
-      std::vector<std::string> iwords = Tokenizer(istr, ",").as_vector();
-      infield = iwords.size();
-
-      jn = strlen(arg[iarg + 2]) + 1;
-      jstr = new char[jn];
-      strcpy(jstr, arg[iarg + 2]);
-      std::vector<std::string> jwords = Tokenizer(jstr, ",").as_vector();
-      jnfield = jwords.size();
-
-      for (i = 0; i < infield; i++) {
-        const char *ifield = iwords[i].c_str();
-        utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error);
-
-        for (j = 0; j < jnfield; j++) {
-          const char *jfield = jwords[j].c_str();
-          utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error);
-
-          for (itype = inlo; itype <= inhi; itype++) {
-            for (jtype = jnlo; jtype <= jnhi; jtype++) {
-              type_filter[itype][jtype] = 1;
-              type_filter[jtype][itype] = 1;
-            }
-          }
-        }
-      }
-
-      delete[] istr;
-      delete[] jstr;
-
-      iarg += 2;
-
-    } else
-      error->all(FLERR, "Invalid keyword in fix pair/tracker command");
-
-    iarg++;
-  }
-
-  if (nvalues == 1)
-    size_local_cols = 0;
-  else
-    size_local_cols = nvalues;
-
-  nmax = 0;
-  ncount = 0;
-  vector = NULL;
-  array = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixPairTracker::~FixPairTracker()
-{
-  delete[] pack_choice;
-
-  memory->destroy(vector);
-  memory->destroy(array);
-  memory->destroy(type_filter);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixPairTracker::setmask()
-{
-  int mask = 0;
-  mask |= POST_FORCE;
-  return mask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::init()
-{
-  // Set size of array/vector
-  ncount = 0;
-
-  if (ncount > nmax) reallocate(ncount);
-
-  size_local_rows = ncount;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::lost_contact(int i, int j, double time_tmp, double nstep_tmp, double rsum_tmp,
-                                  double rmin_tmp)
-{
-
-  double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
-  if ((time - time_tmp) < tmin) return;
-
-  if (type_filter) {
-    int *type = atom->type;
-    if (type_filter[type[i]][type[j]] == 0) return;
-  }
-
-  int *mask = atom->mask;
-  if (!(mask[i] & groupbit)) return;
-  if (!(mask[j] & groupbit)) return;
-
-  if (ncount == nmax) reallocate(ncount);
-
-  index_i = i;
-  index_j = j;
-
-  rmin = rmin_tmp;
-  rsum = rsum_tmp;
-  time_initial = time_tmp;
-  nstep_initial = nstep_tmp;
-
-  // fill vector or array with local values
-  if (nvalues == 1) {
-    (this->*pack_choice[0])(0);
-  } else {
-    for (int k = 0; k < nvalues; k++) { (this->*pack_choice[k])(k); }
-  }
-
-  ncount += 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::post_force(int /*vflag*/)
-{
-  if (update->ntimestep % nevery == 0) {
-    size_local_rows = ncount;
-    ncount = 0;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::reallocate(int n)
-{
-  // grow vector or array
-  while (nmax <= n) nmax += DELTA;
-
-  if (nvalues == 1) {
-    memory->grow(vector, nmax, "fix_pair_tracker:vector");
-    vector_local = vector;
-  } else {
-    memory->grow(array, nmax, nvalues, "fix_pair_tracker:array");
-    array_local = array;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   memory usage of local data
-------------------------------------------------------------------------- */
-
-double FixPairTracker::memory_usage()
-{
-  double bytes = nmax * (double) nvalues * sizeof(double);
-  bytes += nmax * 2 * sizeof(int);
-  return bytes;
-}
-
-/* ----------------------------------------------------------------------
-   one method for every keyword fix pair/tracker can output
-   the atom property is packed into a local vector or array
-------------------------------------------------------------------------- */
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_time_created(int n)
-{
-  if (nvalues == 1)
-    vector[ncount] = time_initial;
-  else
-    array[ncount][n] = time_initial;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_time_broken(int n)
-{
-  double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
-  if (nvalues == 1)
-    vector[ncount] = time;
-  else
-    array[ncount][n] = time;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_time_total(int n)
-{
-  double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
-  if (nvalues == 1)
-    vector[ncount] = time - time_initial;
-  else
-    array[ncount][n] = time - time_initial;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_id1(int n)
-{
-  tagint *tag = atom->tag;
-
-  if (nvalues == 1)
-    vector[ncount] = tag[index_i];
-  else
-    array[ncount][n] = tag[index_i];
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_id2(int n)
-{
-  tagint *tag = atom->tag;
-
-  if (nvalues == 1)
-    vector[ncount] = tag[index_j];
-  else
-    array[ncount][n] = tag[index_j];
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_x(int n)
-{
-  double **x = atom->x;
-
-  if (nvalues == 1)
-    vector[ncount] = (x[index_i][0] + x[index_j][0]) / 2;
-  else
-    array[ncount][n] = (x[index_i][0] + x[index_j][0]) / 2;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_y(int n)
-{
-  double **x = atom->x;
-
-  if (nvalues == 1)
-    vector[ncount] = (x[index_i][1] + x[index_j][1]) / 2;
-  else
-    array[ncount][n] = (x[index_i][1] + x[index_j][1]) / 2;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_z(int n)
-{
-  double **x = atom->x;
-
-  if (nvalues == 1)
-    vector[ncount] = (x[index_i][2] + x[index_j][2]) / 2;
-  else
-    array[ncount][n] = (x[index_i][2] + x[index_j][2]) / 2;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_rmin(int n)
-{
-  if (nvalues == 1)
-    vector[ncount] = rmin;
-  else
-    array[ncount][n] = rmin;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPairTracker::pack_rave(int n)
-{
-  if (nvalues == 1)
-    vector[ncount] = rsum / (update->ntimestep - nstep_initial);
-  else
-    array[ncount][n] = rsum / (update->ntimestep - nstep_initial);
-}

src/MISC/pair_tracker.cpp

@@ -19,20 +19,23 @@
 #include "fix.h"
 #include "fix_dummy.h"
 #include "fix_neigh_history.h"
-#include "fix_pair_tracker.h"
+#include "fix_store_local.h"
 #include "force.h"
 #include "memory.h"
 #include "modify.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "neighbor.h"
+#include "tokenizer.h"
+#include "utils.h"
 #include "update.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
+PairTracker::PairTracker(LAMMPS *lmp) :
+    Pair(lmp), pack_choice(NULL)
 {
   single_enable = 1;
   no_virial_fdotr_compute = 1;
@@ -43,6 +46,7 @@ PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
   nondefault_history_transfer = 1;
 
   finitecutflag = 0;
+  tmin = -1;
 
   // create dummy fix as placeholder for FixNeighHistory
   // this is so final order of Modify:fix will conform to input script
@@ -50,6 +54,12 @@ PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
   fix_history = nullptr;
   modify->add_fix("NEIGH_HISTORY_TRACK_DUMMY all DUMMY");
   fix_dummy = (FixDummy *) modify->fix[modify->nfix - 1];
+
+  fix_store_local = nullptr;
+
+  output_data = nullptr;
+  pack_choice = nullptr;
+  type_filter = nullptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -71,6 +81,12 @@ PairTracker::~PairTracker()
     delete[] maxrad_dynamic;
     delete[] maxrad_frozen;
   }
+
+  delete[] pack_choice;
+
+  delete [] id_fix_store_local;
+  memory->destroy(output_data);
+  memory->destroy(type_filter);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -130,8 +146,9 @@ void PairTracker::compute(int eflag, int vflag)
 
           data = &alldata[size_history * jj];
           if (touch[jj] == 1) {
-            fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]);
+            process_data(i, j, data);
           }
+
           touch[jj] = 0;
           data[0] = 0.0;    // initial time
           data[1] = 0.0;    // initial timestep
@@ -159,7 +176,7 @@ void PairTracker::compute(int eflag, int vflag)
 
           data = &alldata[size_history * jj];
           if (touch[jj] == 1) {
-            fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]);
+            process_data(i, j, data);
           }
 
           touch[jj] = 0;
@@ -216,14 +233,98 @@ void PairTracker::allocate()
 
 void PairTracker::settings(int narg, char **arg)
 {
-  if (narg != 0 && narg != 1) error->all(FLERR, "Illegal pair_style command");
+  if (narg < 1) error->all(FLERR, "Illegal pair_style command");
 
-  if (narg == 1) {
-    if (strcmp(arg[0], "finite") == 0)
+  id_fix_store_local = utils::strdup(arg[0]);
+
+  // If optional arguments included, this will be oversized
+  pack_choice = new FnPtrPack[narg - 1];
+
+  nvalues = 0;
+  int iarg = 1;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "finite") == 0) {
       finitecutflag = 1;
-    else
+    } else if (strcmp(arg[iarg], "id1") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_id1;
+    } else if (strcmp(arg[iarg], "id2") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_id2;
+    } else if (strcmp(arg[iarg], "time/created") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_time_created;
+    } else if (strcmp(arg[iarg], "time/broken") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_time_broken;
+    } else if (strcmp(arg[iarg], "time/total") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_time_total;
+    } else if (strcmp(arg[iarg], "x") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_x;
+    } else if (strcmp(arg[iarg], "y") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_y;
+    } else if (strcmp(arg[iarg], "z") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_z;
+    } else if (strcmp(arg[iarg], "r/min") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_rmin;
+    } else if (strcmp(arg[iarg], "r/ave") == 0) {
+      pack_choice[nvalues++] = &PairTracker::pack_rave;
+    } else if (strcmp(arg[iarg], "time/min") == 0) {
+      if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in pair tracker command");
+      tmin = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      iarg++;
+
+    } else if (strcmp(arg[iarg], "type/include") == 0) {
+      if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in pair tracker command");
+      int ntypes = atom->ntypes;
+
+      int i, j, itype, jtype, in, jn, infield, jnfield;
+      int inlo, inhi, jnlo, jnhi;
+      char *istr, *jstr;
+      if (!type_filter) {
+        memory->create(type_filter, ntypes + 1, ntypes + 1, "pair/tracker:type_filter");
+
+        for (i = 0; i <= ntypes; i++) {
+          for (j = 0; j <= ntypes; j++) type_filter[i][j] = 0;
+        }
+      }
+
+      in = strlen(arg[iarg + 1]) + 1;
+      istr = new char[in];
+      strcpy(istr, arg[iarg + 1]);
+      std::vector<std::string> iwords = Tokenizer(istr, ",").as_vector();
+      infield = iwords.size();
+
+      jn = strlen(arg[iarg + 2]) + 1;
+      jstr = new char[jn];
+      strcpy(jstr, arg[iarg + 2]);
+      std::vector<std::string> jwords = Tokenizer(jstr, ",").as_vector();
+      jnfield = jwords.size();
+
+      for (i = 0; i < infield; i++) {
+        const char *ifield = iwords[i].c_str();
+        utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error);
+
+        for (j = 0; j < jnfield; j++) {
+          const char *jfield = jwords[j].c_str();
+          utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error);
+
+          for (itype = inlo; itype <= inhi; itype++) {
+            for (jtype = jnlo; jtype <= jnhi; jtype++) {
+              type_filter[itype][jtype] = 1;
+              type_filter[jtype][itype] = 1;
+            }
+          }
+        }
+      }
+      delete[] istr;
+      delete[] jstr;
+      iarg += 2;
+
+    } else {
       error->all(FLERR, "Illegal pair_style command");
+    }
+    iarg ++;
   }
+
+  if (nvalues == 0) error->all(FLERR, "Must request at least one value to output");
+  memory->create(output_data, nvalues, "pair/tracker:output_data");
 }
 
 /* ----------------------------------------------------------------------
@@ -262,7 +363,6 @@ void PairTracker::coeff(int narg, char **arg)
 void PairTracker::init_style()
 {
   int i;
-
   // error and warning checks
 
   if (!atom->radius_flag && finitecutflag)
@@ -275,7 +375,7 @@ void PairTracker::init_style()
     neighbor->requests[irequest]->size = 1;
     neighbor->requests[irequest]->history = 1;
     // history flag won't affect results, but match granular pairstyles
-    // so neighborlist can be copied to reduce overhead
+    // to copy neighbor list and reduce overhead
   }
 
   // if history is stored and first init, create Fix to store history
@@ -284,7 +384,7 @@ void PairTracker::init_style()
 
   if (fix_history == nullptr) {
     modify->replace_fix("NEIGH_HISTORY_TRACK_DUMMY",
-                        fmt::format("NEIGH_HISTORY_TRACK all NEIGH_HISTORY {}", size_history), 1);
+        fmt::format("NEIGH_HISTORY_TRACK all NEIGH_HISTORY {}", size_history), 1);
     int ifix = modify->find_fix("NEIGH_HISTORY_TRACK");
     fix_history = (FixNeighHistory *) modify->fix[ifix];
     fix_history->pair = this;
@@ -295,9 +395,14 @@ void PairTracker::init_style()
 
     if (force->pair->beyond_contact)
       error->all(FLERR,
-                 "Pair tracker incompatible with granular pairstyles that extend beyond contact");
-    // check for FixPour and FixDeposit so can extract particle radii
+        "Pair tracker incompatible with granular pairstyles that extend beyond contact");
 
+    // check for FixFreeze and set freeze_group_bit
+    int ifix = modify->find_fix_by_style("^freeze");
+    if (ifix < 0) freeze_group_bit = 0;
+    else freeze_group_bit = modify->fix[ifix]->groupbit;
+
+    // check for FixPour and FixDeposit so can extract particle radii
     int ipour;
     for (ipour = 0; ipour < modify->nfix; ipour++)
       if (strcmp(modify->fix[ipour]->style, "pour") == 0) break;
@@ -310,7 +415,6 @@ void PairTracker::init_style()
 
     // set maxrad_dynamic and maxrad_frozen for each type
     // include future FixPour and FixDeposit particles as dynamic
-
     int itype;
     for (i = 1; i <= atom->ntypes; i++) {
       onerad_dynamic[i] = onerad_frozen[i] = 0.0;
@@ -343,9 +447,13 @@ void PairTracker::init_style()
   if (ifix < 0) error->all(FLERR, "Could not find pair fix neigh history ID");
   fix_history = (FixNeighHistory *) modify->fix[ifix];
 
-  ifix = modify->find_fix_by_style("pair/tracker");
-  if (ifix < 0) error->all(FLERR, "Cannot use pair tracker without fix pair/tracker");
-  fix_pair_tracker = (FixPairTracker *) modify->fix[ifix];
+  ifix = modify->find_fix(id_fix_store_local);
+  if (ifix < 0) error->all(FLERR, "Cannot find fix store/local");
+  if (strcmp(modify->fix[ifix]->style, "store/local") != 0)
+      error->all(FLERR, "Incorrect fix style matched, not store/local");
+  fix_store_local = (FixStoreLocal *) modify->fix[ifix];
+  if (fix_store_local->nvalues != nvalues)
+    error->all(FLERR, "Inconsistent number of output variables in fix store/local");
 }
 
 /* ----------------------------------------------------------------------
@@ -357,7 +465,6 @@ double PairTracker::init_one(int i, int j)
   if (!allocated) allocate();
 
   // always mix prefactors geometrically
-
   if (setflag[i][j] == 0) { cut[i][j] = mix_distance(cut[i][i], cut[j][j]); }
 
   cut[j][i] = cut[i][j];
@@ -469,3 +576,105 @@ double PairTracker::radii2cut(double r1, double r2)
   double cut = r1 + r2;
   return cut;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::process_data(int i, int j, double * input_data)
+{
+  double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
+  int time_initial = (int) input_data[0];
+  if ((time - time_initial) < tmin) return;
+
+  if (type_filter) {
+    int *type = atom->type;
+    if (type_filter[type[i]][type[j]] == 0) return;
+  }
+
+  for (int k = 0; k < nvalues; k++) (this->*pack_choice[k])(k, i, j, input_data);
+  fix_store_local->add_data(output_data, i, j);
+}
+
+/* ----------------------------------------------------------------------
+   one method for every keyword fix pair/tracker can output
+   the atom property is packed into a local vector or array
+------------------------------------------------------------------------- */
+
+void PairTracker::pack_time_created(int n, int i, int j, double * data)
+{
+  double time_initial = data[0];
+  output_data[n] = time_initial;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_time_broken(int n, int i, int j, double * data)
+{
+  double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
+  output_data[n] = time;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_time_total(int n, int i, int j, double * data)
+{
+  double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
+  double time_initial = data[0];  
+  output_data[n] = time - time_initial;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_id1(int n, int i, int j, double * data)
+{
+  tagint *tag = atom->tag;
+  output_data[n] = tag[i];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_id2(int n, int i, int j, double * data)
+{
+  tagint *tag = atom->tag;
+  output_data[n] = tag[j];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_x(int n, int i, int j, double * data)
+{
+  double **x = atom->x;
+  output_data[n] = (x[i][0] + x[j][0])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_y(int n, int i, int j, double * data)
+{
+  double **x = atom->x;
+  output_data[n] = (x[i][1] + x[j][1])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_z(int n, int i, int j, double * data)
+{
+  double **x = atom->x;
+  output_data[n] = (x[i][2] + x[j][2])*0.5;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_rmin(int n, int i, int j, double * data)
+{
+  double rmin = data[3];
+  output_data[n] = rmin;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::pack_rave(int n, int i, int j, double * data)
+{
+  double rsum = data[2];
+  double nstep_initial = data[1];
+  output_data[n] = rsum / (update->ntimestep - nstep_initial);
+}

src/MISC/pair_tracker.h

@@ -51,11 +51,33 @@ class PairTracker : public Pair {
   double *maxrad_dynamic, *maxrad_frozen;
   int freeze_group_bit;
 
+  char *id_fix_store_local;
   class FixDummy *fix_dummy;
   class FixNeighHistory *fix_history;
-  class FixPairTracker *fix_pair_tracker;
+  class FixStoreLocal *fix_store_local;
 
+  int **type_filter;
+  double tmin;
+
+  int nvalues, ncount;
+  double *output_data;
+  typedef void (PairTracker::*FnPtrPack)(int, int, int, double *);
+  FnPtrPack *pack_choice;    // ptrs to pack functions
+
+  void pack_id1(int, int, int, double *);
+  void pack_id2(int, int, int, double *);
+  void pack_time_created(int, int, int, double *);
+  void pack_time_broken(int, int, int, double *);
+  void pack_time_total(int, int, int, double *);
+  void pack_x(int, int, int, double *);
+  void pack_y(int, int, int, double *);
+  void pack_z(int, int, int, double *);
+  void pack_rmin(int, int, int, double *);
+  void pack_rave(int, int, int, double *);
+
+  void transfer_data(int, int);
   void transfer_history(double *, double *);
+  void process_data(int, int, double *);
   void allocate();
 };
 
@@ -72,20 +94,37 @@ Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
+E: Invalid keyword in pair tracker command
+
+Self-explanatory.
+
 E: Incorrect args for pair coefficients
 
 Self-explanatory.  Check the input script or data file.
 
+E: Must request at least one value to output
+
+Must include at least one bond property to store in fix store/local
+
 E: Pair tracker requires atom attribute radius for finite cutoffs
 
 The atom style defined does not have these attributes.
 
+E: Pair tracker incompatible with granular pairstyles that extend beyond contact
+
+Self-explanatory.
+
 E: Could not find pair fix neigh history ID
 
 The associated fix neigh/history is missing
 
-E: Cannot use pair tracker without fix pair/tracker
+E: Cannot use pair tracker without fix store/local
 
-This pairstyle requires one to define a pair/tracker fix
+The associated fix store/local does not exist
+
+E: Inconsistent number of output variables in fix store/local
+
+The number of values specified in fix store/local disagrees with
+the number of values requested in pair tracker
 
 */

src/MOLECULE/bond_quartic.cpp

@@ -34,6 +34,7 @@ using namespace LAMMPS_NS;
 
 BondQuartic::BondQuartic(LAMMPS *lmp) : Bond(lmp)
 {
+  partial_flag = 1;
   TWO_1_3 = pow(2.0,(1.0/3.0));
 }
 

src/Makefile

@@ -53,6 +53,7 @@ PACKAGE = \
 	awpmd \
 	bocs \
 	body \
+    bpm \
 	brownian \
 	cg-dna \
 	cg-sdk \
@@ -147,7 +148,8 @@ PACKMOST = \
 	asphere \
 	bocs \
 	body \
-	brownian \
+	bpm
+    brownian \
 	cg-dna \
 	cg-sdk \
 	class2 \

src/OPENMP/fix_neigh_history_omp.cpp

@@ -576,7 +576,7 @@ void FixNeighHistoryOMP::post_neighbor()
 
       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
-        rflag = sbmask(j);
+        rflag = histmask(j);
         j &= NEIGHMASK;
         jlist[jj] = j;
 

src/OPENMP/npair_half_size_bin_newtoff_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_bin_newtoff_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 #include "npair_omp.h"
@@ -40,8 +43,10 @@ NPairHalfSizeBinNewtoffOmp::NPairHalfSizeBinNewtoffOmp(LAMMPS *lmp) :
 void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 
@@ -50,7 +55,8 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,ibin;
+  int i,j,jh,k,n,ibin,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -61,7 +67,14 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -81,6 +94,11 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
     ztmp = x[i][2];
     radi = radius[i];
     ibin = atom2bin[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in surrounding bins in stencil including self
     // only store pair if i < j
@@ -100,10 +118,23 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >=0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/OPENMP/npair_half_size_bin_newton_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_bin_newton_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 #include "npair_omp.h"
@@ -39,8 +42,10 @@ NPairHalfSizeBinNewtonOmp::NPairHalfSizeBinNewtonOmp(LAMMPS *lmp) :
 void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 
@@ -49,7 +54,8 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,ibin;
+  int i,j,jh,k,n,ibin,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -58,7 +64,14 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -77,6 +90,11 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over rest of atoms in i's bin, ghosts are at end of linked list
     // if j is owned atom, store it, since j is beyond i in linked list
@@ -101,10 +119,23 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >=0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 
@@ -123,10 +154,23 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >=0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_bin_newton_tri_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 #include "npair_omp.h"
@@ -39,8 +42,10 @@ NPairHalfSizeBinNewtonTriOmp::NPairHalfSizeBinNewtonTriOmp(LAMMPS *lmp) :
 void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -48,7 +53,8 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,ibin;
+  int i,j,jh,k,n,ibin,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -59,7 +65,14 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -78,6 +91,11 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in bins in stencil
     // pairs for atoms j "below" i are excluded
@@ -107,10 +125,23 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >=0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/OPENMP/npair_half_size_multi_newtoff_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_multi_newtoff_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -41,8 +44,10 @@ NPairHalfSizeMultiNewtoffOmp::NPairHalfSizeMultiNewtoffOmp(LAMMPS *lmp) : NPair(
 void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -50,7 +55,9 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -63,7 +70,14 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -84,6 +98,11 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     ibin = atom2bin[i];
 
@@ -104,27 +123,40 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
       ns = nstencil_multi[icollection][jcollection];
 
       for (k = 0; k < ns; k++) {
-            js = binhead_multi[jcollection][jbin + s[k]];
-            for (j = js; j >=0; j = bins[j]) {
-              if (j <= i) continue;
+        js = binhead_multi[jcollection][jbin + s[k]];
+        for (j = js; j >=0; j = bins[j]) {
+          if (j <= i) continue;
 
           jtype = type[j];
-              if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
-              delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                  neighptr[n++] = j ^ mask_history;
-                else
-                  neighptr[n++] = j;
-              }
-            }
+          if (rsq <= cutdistsq) {
+            jh = j;
+            if (history && rsq < radsum*radsum)
+              jh = jh ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >=0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = jh;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = jh;
+              else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            } else neighptr[n++] = jh;
+          }
+        }
       }
     }
 

src/OPENMP/npair_half_size_multi_newton_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_multi_newton_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -40,8 +43,10 @@ NPairHalfSizeMultiNewtonOmp::NPairHalfSizeMultiNewtonOmp(LAMMPS *lmp) : NPair(lm
 void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -49,7 +54,9 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -62,7 +69,14 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -83,6 +97,11 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     ibin = atom2bin[i];
 
@@ -107,7 +126,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
         //   if j is owned atom, store it if j > i
         //   if j is ghost, only store if j coords are "above and to the right" of i
 
-            for (j = js; j >= 0; j = bins[j]) {
+        for (j = js; j >= 0; j = bins[j]) {
           if(icollection != jcollection and j < i) continue;
 
               if (j >= nlocal) {
@@ -121,19 +140,32 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
           jtype = type[j];
           if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
-              delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                  neighptr[n++] = j ^ mask_history;
-                else
-                  neighptr[n++] = j;
-              }
+          if (rsq <= cutdistsq) {
+            jh = j;
+            if (history && rsq < radsum*radsum)
+              jh = jh ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >=0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = jh;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = jh;
+              else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            } else neighptr[n++] = jh;
+          }
         }
       }
 
@@ -142,31 +174,43 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
       // stencil is half if i same size as j
       // stencil is full if i smaller than j
 
-          s = stencil_multi[icollection][jcollection];
-          ns = nstencil_multi[icollection][jcollection];
+      s = stencil_multi[icollection][jcollection];
+      ns = nstencil_multi[icollection][jcollection];
 
-          for (k = 0; k < ns; k++) {
-            js = binhead_multi[jcollection][jbin + s[k]];
-            for (j = js; j >= 0; j = bins[j]) {
+      for (k = 0; k < ns; k++) {
+        js = binhead_multi[jcollection][jbin + s[k]];
+        for (j = js; j >= 0; j = bins[j]) {
 
           jtype = type[j];
-              if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
           delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                    neighptr[n++] = j ^ mask_history;
-                else
-                    neighptr[n++] = j;
-              }
-            }
+          if (rsq <= cutdistsq) {
+            if (history && rsq < radsum*radsum)
+              j = j ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >=0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = j;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = j;
+              else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+            } else neighptr[n++] = j;
           }
+        }
+      }
     }
 
     ilist[i] = i;

src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_multi_newton_tri_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -41,8 +44,10 @@ NPairHalfSizeMultiNewtonTriOmp::NPairHalfSizeMultiNewtonTriOmp(LAMMPS *lmp) :
 void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -50,7 +55,9 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -63,7 +70,14 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -84,6 +98,11 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     ibin = atom2bin[i];
 
@@ -104,12 +123,12 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
       //         (equal zyx and j <= i)
       // latter excludes self-self interaction but allows superposed atoms
 
-          s = stencil_multi[icollection][jcollection];
-          ns = nstencil_multi[icollection][jcollection];
+      s = stencil_multi[icollection][jcollection];
+      ns = nstencil_multi[icollection][jcollection];
 
-          for (k = 0; k < ns; k++) {
-            js = binhead_multi[jcollection][jbin + s[k]];
-            for (j = js; j >= 0; j = bins[j]) {
+      for (k = 0; k < ns; k++) {
+        js = binhead_multi[jcollection][jbin + s[k]];
+        for (j = js; j >= 0; j = bins[j]) {
 
           // if same size (same collection), use half stencil
           if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@@ -126,21 +145,34 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
           jtype = type[j];
           if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
-              delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                    neighptr[n++] = j ^ mask_history;
-                else
-                    neighptr[n++] = j;
-              }
-            }
+          if (rsq <= cutdistsq) {
+            jh = j;
+            if (history && rsq < radsum*radsum)
+              jh = jh ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >=0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = jh;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = jh;
+              else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            } else neighptr[n++] = jh;
           }
+        }
+      }
     }
 
     ilist[i] = i;

src/OPENMP/npair_half_size_multi_old_newtoff_omp.cpp

@@ -18,6 +18,8 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
+#include "domain.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -40,8 +42,10 @@ NPairHalfSizeMultiOldNewtoffOmp::NPairHalfSizeMultiOldNewtoffOmp(LAMMPS *lmp) :
 void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -49,7 +53,8 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,itype,jtype,ibin,ns;
+  int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -59,7 +64,14 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -79,6 +91,11 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in other bins in stencil including self
     // only store pair if i < j
@@ -107,10 +124,23 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list)
         cutdistsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutdistsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >=0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/OPENMP/npair_half_size_multi_old_newton_omp.cpp

@@ -18,6 +18,8 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
+#include "domain.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -39,8 +41,10 @@ NPairHalfSizeMultiOldNewtonOmp::NPairHalfSizeMultiOldNewtonOmp(LAMMPS *lmp) :
 void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -48,7 +52,8 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,itype,jtype,ibin,ns;
+  int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -58,7 +63,14 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -78,6 +90,11 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over rest of atoms in i's bin, ghosts are at end of linked list
     // if j is owned atom, store it, since j is beyond i in linked list
@@ -102,10 +119,23 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list)
       cutdistsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutdistsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >=0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 
@@ -132,10 +162,23 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list)
         cutdistsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutdistsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >=0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp

@@ -18,6 +18,8 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
+#include "domain.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -39,8 +41,10 @@ NPairHalfSizeMultiOldNewtonTriOmp::NPairHalfSizeMultiOldNewtonTriOmp(LAMMPS *lmp
 void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 #if defined(_OPENMP)
@@ -48,7 +52,8 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,k,n,itype,jtype,ibin,ns;
+  int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -58,7 +63,14 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -78,6 +90,11 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in bins, including self, in stencil
     // skip if i,j neighbor cutoff is less than bin distance
@@ -115,10 +132,23 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
         cutdistsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutdistsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >=0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/OPENMP/npair_half_size_nsq_newtoff_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_nsq_newtoff_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "group.h"
 #include "my_page.h"
 #include "neigh_list.h"
@@ -42,8 +45,10 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
   const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 
@@ -52,7 +57,8 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,n;
+  int i,j,jh,n,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -61,9 +67,17 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
   int nall = atom->nlocal + atom->nghost;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
@@ -81,6 +95,11 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over remaining atoms, owned and ghost
 
@@ -96,10 +115,23 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >=0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 

src/OPENMP/npair_half_size_nsq_newton_omp.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_nsq_newton_omp.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "group.h"
 #include "my_page.h"
 #include "neigh_list.h"
@@ -42,9 +45,11 @@ NPairHalfSizeNsqNewtonOmp::NPairHalfSizeNsqNewtonOmp(LAMMPS *lmp) :
 void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
 {
   const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
-  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;;
+  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
+  const int molecular = atom->molecular;
+  const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
   const int history = list->history;
-  const int mask_history = 3 << SBBITS;
+  const int mask_history = 1 << HISTBITS;
 
   NPAIR_OMP_INIT;
 
@@ -53,7 +58,8 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
 #endif
   NPAIR_OMP_SETUP(nlocal);
 
-  int i,j,n,itag,jtag;
+  int i,j,jh,n,itag,jtag,which,imol,iatom;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -64,6 +70,13 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
   int *type = atom->type;
   int *mask = atom->mask;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
+
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+
   int nall = atom->nlocal + atom->nghost;
 
   int *ilist = list->ilist;
@@ -84,6 +97,11 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over remaining atoms, owned and ghost
 
@@ -115,10 +133,23 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >=0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 

src/atom.cpp

@@ -122,6 +122,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   omega = angmom = torque = nullptr;
   radius = rmass = nullptr;
   ellipsoid = line = tri = body = nullptr;
+  quat = nullptr;
 
   // molecular systems
 
@@ -419,6 +420,8 @@ void Atom::peratom_create()
   add_peratom("tri",&tri,INT,0);
   add_peratom("body",&body,INT,0);
 
+  add_peratom("quat",&quat,DOUBLE,4);
+
   // MOLECULE package
 
   add_peratom("molecule",&molecule,tagintsize,0);
@@ -639,6 +642,7 @@ void Atom::set_atomflag_defaults()
   // identical list as 2nd customization in atom.h
 
   sphere_flag = ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
+  quat_flag = 0;
   peri_flag = electron_flag = 0;
   wavepacket_flag = sph_flag = 0;
   molecule_flag = molindex_flag = molatom_flag = 0;
@@ -2660,6 +2664,10 @@ length of the data area, and a short description.
      - int
      - 1
      - 1 if the particle is a body particle, 0 if not
+  * - quat
+     - double
+     - 4
+     - four quaternion components of the particles
    * - i_name
      - int
      - 1
@@ -2715,6 +2723,7 @@ void *Atom::extract(const char *name)
   if (strcmp(name,"line") == 0) return (void *) line;
   if (strcmp(name,"tri") == 0) return (void *) tri;
   if (strcmp(name,"body") == 0) return (void *) body;
+  if (strcmp(name,"quat") == 0) return (void *) quat;
 
   if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
   if (strcmp(name,"s0") == 0) return (void *) s0;
@@ -2837,6 +2846,7 @@ int Atom::extract_datatype(const char *name)
   if (strcmp(name,"line") == 0) return LAMMPS_INT;
   if (strcmp(name,"tri") == 0) return LAMMPS_INT;
   if (strcmp(name,"body") == 0) return LAMMPS_INT;
+  if (strcmp(name,"quat") == 0) return LAMMPS_DOUBLE_2D;
 
   if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;

src/atom.h

@@ -79,6 +79,7 @@ class Atom : protected Pointers {
   double *radius;
   double **omega, **angmom, **torque;
   int *ellipsoid, *line, *tri, *body;
+  double **quat;
 
   // molecular systems
 
@@ -180,7 +181,7 @@ class Atom : protected Pointers {
 
   int molecule_flag, molindex_flag, molatom_flag;
   int q_flag, mu_flag;
-  int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag;
+  int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag, quat_flag;
   int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag;
   int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag;
   int rho_flag, esph_flag, cv_flag, vest_flag;

src/atom_vec_sphere.cpp

@@ -6,7 +6,7 @@
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributead under
+   certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.

src/bond.cpp

@@ -42,8 +42,11 @@ Bond::Bond(LAMMPS *lmp) : Pointers(lmp)
   virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
   writedata = 1;
 
+  comm_forward = comm_reverse = comm_reverse_off = 0;
+
   allocated = 0;
   suffix_flag = Suffix::NONE;
+  partial_flag = 0;
 
   maxeatom = maxvatom = 0;
   eatom = nullptr;

src/bond.h

@@ -25,11 +25,16 @@ class Bond : protected Pointers {
  public:
   int allocated;
   int *setflag;
+  int partial_flag;          // 1 if bond type can be set to 0 and deleted
   int writedata;             // 1 if writes coeffs to data file
   double energy;             // accumulated energies
   double virial[6];          // accumulated virial: xx,yy,zz,xy,xz,yz
   double *eatom, **vatom;    // accumulated per-atom energy/virial
 
+  int comm_forward;          // size of forward communication (0 if none)
+  int comm_reverse;          // size of reverse communication (0 if none)
+  int comm_reverse_off;      // size of reverse comm even if newton off
+
   int reinitflag;    // 1 if compatible with fix adapt and alike
 
   // KOKKOS host/device flag and data masks
@@ -55,6 +60,10 @@ class Bond : protected Pointers {
   virtual double memory_usage();
   virtual void *extract(const char *, int &) { return nullptr; }
   virtual void reinit();
+  virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;}
+  virtual void unpack_forward_comm(int, int, double *) {}
+  virtual int pack_reverse_comm(int, int, double *) {return 0;}
+  virtual void unpack_reverse_comm(int, int *, double *) {}
 
   void write_file(int, char **);
 

src/comm.cpp

@@ -217,6 +217,9 @@ void Comm::init()
   if (force->pair) maxforward = MAX(maxforward,force->pair->comm_forward);
   if (force->pair) maxreverse = MAX(maxreverse,force->pair->comm_reverse);
 
+  if (force->bond) maxforward = MAX(maxforward,force->bond->comm_forward);
+  if (force->bond) maxreverse = MAX(maxreverse,force->bond->comm_reverse);
+
   for (int i = 0; i < modify->nfix; i++) {
     maxforward = MAX(maxforward,modify->fix[i]->comm_forward);
     maxreverse = MAX(maxreverse,modify->fix[i]->comm_reverse);
@@ -234,6 +237,7 @@ void Comm::init()
 
   if (force->newton == 0) maxreverse = 0;
   if (force->pair) maxreverse = MAX(maxreverse,force->pair->comm_reverse_off);
+  if (force->bond) maxreverse = MAX(maxreverse,force->bond->comm_reverse_off);
 
   // maxexchange_atom = size of an exchanged atom, set by AtomVec
   //   only needs to be set if size > BUFEXTRA

src/comm.h

@@ -84,6 +84,8 @@ class Comm : protected Pointers {
 
   virtual void forward_comm_pair(class Pair *) = 0;
   virtual void reverse_comm_pair(class Pair *) = 0;
+  virtual void forward_comm_bond(class Bond *) = 0;
+  virtual void reverse_comm_bond(class Bond *) = 0;
   virtual void forward_comm_fix(class Fix *, int size = 0) = 0;
   virtual void reverse_comm_fix(class Fix *, int size = 0) = 0;
   virtual void reverse_comm_fix_variable(class Fix *) = 0;

src/comm_brick.cpp

@@ -20,6 +20,7 @@
 
 #include "atom.h"
 #include "atom_vec.h"
+#include "bond.h"
 #include "compute.h"
 #include "domain.h"
 #include "dump.h"
@@ -1074,6 +1075,79 @@ void CommBrick::reverse_comm_pair(Pair *pair)
   }
 }
 
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Bond
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommBrick::forward_comm_bond(Bond *bond)
+{
+  int iswap,n;
+  double *buf;
+  MPI_Request request;
+
+  int nsize = bond->comm_forward;
+
+  for (iswap = 0; iswap < nswap; iswap++) {
+
+    // pack buffer
+
+    n = bond->pack_forward_comm(sendnum[iswap],sendlist[iswap],
+                                buf_send,pbc_flag[iswap],pbc[iswap]);
+
+    // exchange with another proc
+    // if self, set recv buffer to send buffer
+
+    if (sendproc[iswap] != me) {
+      if (recvnum[iswap])
+        MPI_Irecv(buf_recv,nsize*recvnum[iswap],MPI_DOUBLE,
+                  recvproc[iswap],0,world,&request);
+      if (sendnum[iswap])
+        MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
+      if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
+      buf = buf_recv;
+    } else buf = buf_send;
+
+    // unpack buffer
+
+    bond->unpack_forward_comm(recvnum[iswap],firstrecv[iswap],buf);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Bond
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommBrick::reverse_comm_bond(Bond *bond)
+{
+  int iswap,n;
+  double *buf;
+  MPI_Request request;
+
+  int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off);
+
+  for (iswap = nswap-1; iswap >= 0; iswap--) {
+
+    // pack buffer
+
+    n = bond->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);
+
+    // exchange with another proc
+    // if self, set recv buffer to send buffer
+
+    if (sendproc[iswap] != me) {
+      if (sendnum[iswap])
+        MPI_Irecv(buf_recv,nsize*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,
+                  world,&request);
+      if (recvnum[iswap])
+        MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world);
+      if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
+      buf = buf_recv;
+    } else buf = buf_send;
+  }
+}
+
 /* ----------------------------------------------------------------------
    forward communication invoked by a Fix
    size/nsize used only to set recv buffer limit

src/comm_brick.h

@@ -33,6 +33,8 @@ class CommBrick : public Comm {
 
   virtual void forward_comm_pair(class Pair *);    // forward comm from a Pair
   virtual void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
+  virtual void forward_comm_bond(class Bond *);    // forward comm from a Bond
+  virtual void reverse_comm_bond(class Bond *);    // reverse comm from a Bond
   virtual void forward_comm_fix(class Fix *, int size = 0);
   // forward comm from a Fix
   virtual void reverse_comm_fix(class Fix *, int size = 0);

src/comm_tiled.cpp

@@ -21,6 +21,7 @@
 
 #include "atom.h"
 #include "atom_vec.h"
+#include "bond.h"
 #include "compute.h"
 #include "domain.h"
 #include "dump.h"
@@ -1464,6 +1465,99 @@ void CommTiled::reverse_comm_pair(Pair *pair)
   }
 }
 
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Bond
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiled::forward_comm_bond(Bond *bond)
+{
+  int i,irecv,n,nsend,nrecv;
+
+  int nsize = bond->comm_forward;
+
+  for (int iswap = 0; iswap < nswap; iswap++) {
+    nsend = nsendproc[iswap] - sendself[iswap];
+    nrecv = nrecvproc[iswap] - sendself[iswap];
+
+    if (recvother[iswap]) {
+      for (i = 0; i < nrecv; i++)
+        MPI_Irecv(&buf_recv[nsize*forward_recv_offset[iswap][i]],
+                  nsize*recvnum[iswap][i],
+                  MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
+    }
+
+    if (sendother[iswap]) {
+      for (i = 0; i < nsend; i++) {
+        n = bond->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i],
+                                    buf_send,pbc_flag[iswap][i],pbc[iswap][i]);
+        MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][i],0,world);
+      }
+    }
+
+    if (sendself[iswap]) {
+      bond->pack_forward_comm(sendnum[iswap][nsend],sendlist[iswap][nsend],
+                              buf_send,pbc_flag[iswap][nsend],
+                              pbc[iswap][nsend]);
+      bond->unpack_forward_comm(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
+                                buf_send);
+    }
+    if (recvother[iswap]) {
+      for (i = 0; i < nrecv; i++) {
+        MPI_Waitany(nrecv,requests,&irecv,MPI_STATUS_IGNORE);
+        bond->unpack_forward_comm(recvnum[iswap][irecv],firstrecv[iswap][irecv],
+                                  &buf_recv[nsize*
+                                            forward_recv_offset[iswap][irecv]]);
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Bond
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiled::reverse_comm_bond(Bond *bond)
+{
+  int i,irecv,n,nsend,nrecv;
+
+  int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off);
+
+  for (int iswap = nswap-1; iswap >= 0; iswap--) {
+    nsend = nsendproc[iswap] - sendself[iswap];
+    nrecv = nrecvproc[iswap] - sendself[iswap];
+
+    if (sendother[iswap]) {
+      for (i = 0; i < nsend; i++)
+        MPI_Irecv(&buf_recv[nsize*reverse_recv_offset[iswap][i]],
+                  nsize*sendnum[iswap][i],MPI_DOUBLE,
+                  sendproc[iswap][i],0,world,&requests[i]);
+    }
+    if (recvother[iswap]) {
+      for (i = 0; i < nrecv; i++) {
+        n = bond->pack_reverse_comm(recvnum[iswap][i],firstrecv[iswap][i],
+                                    buf_send);
+        MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap][i],0,world);
+      }
+    }
+    if (sendself[iswap]) {
+      bond->pack_reverse_comm(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
+                              buf_send);
+      bond->unpack_reverse_comm(sendnum[iswap][nsend],sendlist[iswap][nsend],
+                                buf_send);
+    }
+    if (sendother[iswap]) {
+      for (i = 0; i < nsend; i++) {
+        MPI_Waitany(nsend,requests,&irecv,MPI_STATUS_IGNORE);
+        bond->unpack_reverse_comm(sendnum[iswap][irecv],sendlist[iswap][irecv],
+                                  &buf_recv[nsize*
+                                            reverse_recv_offset[iswap][irecv]]);
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    forward communication invoked by a Fix
    size/nsize used only to set recv buffer limit

src/comm_tiled.h

@@ -33,6 +33,8 @@ class CommTiled : public Comm {
 
   virtual void forward_comm_pair(class Pair *);    // forward comm from a Pair
   virtual void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
+  virtual void forward_comm_bond(class Bond *);    // forward comm from a Bond
+  virtual void reverse_comm_bond(class Bond *);    // reverse comm from a Bond
   virtual void forward_comm_fix(class Fix *, int size = 0);
   // forward comm from a Fix
   virtual void reverse_comm_fix(class Fix *, int size = 0);

src/compute_property_atom.cpp

@@ -231,25 +231,25 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"quatw") == 0) {
       avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
       avec_body = (AtomVecBody *) atom->style_match("body");
-      if (!avec_ellipsoid && !avec_body)
+      if (!avec_ellipsoid && !avec_body && !atom->quat_flag)
         error->all(FLERR,"Compute property/atom for atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_quatw;
     } else if (strcmp(arg[iarg],"quati") == 0) {
       avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
       avec_body = (AtomVecBody *) atom->style_match("body");
-      if (!avec_ellipsoid && !avec_body)
+      if (!avec_ellipsoid && !avec_body && !atom->quat_flag)
         error->all(FLERR,"Compute property/atom for atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_quati;
     } else if (strcmp(arg[iarg],"quatj") == 0) {
       avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
       avec_body = (AtomVecBody *) atom->style_match("body");
-      if (!avec_ellipsoid && !avec_body)
+      if (!avec_ellipsoid && !avec_body && !atom->quat_flag)
         error->all(FLERR,"Compute property/atom for atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_quatj;
     } else if (strcmp(arg[iarg],"quatk") == 0) {
       avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
       avec_body = (AtomVecBody *) atom->style_match("body");
-      if (!avec_ellipsoid && !avec_body)
+      if (!avec_ellipsoid && !avec_body && !atom->quat_flag)
         error->all(FLERR,"Compute property/atom for atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_quatk;
 
@@ -1334,7 +1334,7 @@ void ComputePropertyAtom::pack_quatw(int n)
       n += nvalues;
     }
 
-  } else {
+  } else if (avec_body) {
     AtomVecBody::Bonus *bonus = avec_body->bonus;
     int *body = atom->body;
     int *mask = atom->mask;
@@ -1346,6 +1346,17 @@ void ComputePropertyAtom::pack_quatw(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+  } else {
+    double **quat = atom->quat;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit)
+        buf[n] = quat[i][0];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
   }
 }
 
@@ -1366,7 +1377,7 @@ void ComputePropertyAtom::pack_quati(int n)
       n += nvalues;
     }
 
-  } else {
+  } else if (avec_body) {
     AtomVecBody::Bonus *bonus = avec_body->bonus;
     int *body = atom->body;
     int *mask = atom->mask;
@@ -1378,6 +1389,17 @@ void ComputePropertyAtom::pack_quati(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+  } else {
+    double **quat = atom->quat;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit)
+        buf[n] = quat[i][1];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
   }
 }
 
@@ -1398,7 +1420,7 @@ void ComputePropertyAtom::pack_quatj(int n)
       n += nvalues;
     }
 
-  } else {
+  } else if (avec_body) {
     AtomVecBody::Bonus *bonus = avec_body->bonus;
     int *body = atom->body;
     int *mask = atom->mask;
@@ -1410,6 +1432,17 @@ void ComputePropertyAtom::pack_quatj(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+  } else {
+    double **quat = atom->quat;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit)
+        buf[n] = quat[i][2];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
   }
 }
 
@@ -1430,7 +1463,7 @@ void ComputePropertyAtom::pack_quatk(int n)
       n += nvalues;
     }
 
-  } else {
+  } else if (avec_body) {
     AtomVecBody::Bonus *bonus = avec_body->bonus;
     int *body = atom->body;
     int *mask = atom->mask;
@@ -1442,6 +1475,17 @@ void ComputePropertyAtom::pack_quatk(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+  } else {
+    double **quat = atom->quat;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit)
+        buf[n] = quat[i][3];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
   }
 }
 

src/fix_neigh_history.cpp

@@ -627,7 +627,7 @@ void FixNeighHistory::post_neighbor()
       j = jlist[jj];
 
       if (use_bit_flag) {
-        rflag = sbmask(j) | pair->beyond_contact;
+        rflag = histmask(j) | pair->beyond_contact;
         j &= NEIGHMASK;
         jlist[jj] = j;
       } else {

src/fix_neigh_history.h

@@ -93,7 +93,8 @@ class FixNeighHistory : public Fix {
   virtual void pre_exchange_no_newton();
   void allocate_pages();
 
-  inline int sbmask(int j) const { return j >> SBBITS & 3; }
+  // Shift by HISTBITS and check the first bit
+  inline int histmask(int j) const { return j >> HISTBITS & 1; }
 };
 
 }    // namespace LAMMPS_NS

src/fix_property_atom.cpp

@@ -550,18 +550,24 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/)
       atom->q[j] = atom->q[i];
     else if (style[nv] == RMASS)
       atom->rmass[j] = atom->rmass[i];
-    else if (style[nv] == IVEC)
+    else if (style[nv] == IVEC) {
       atom->ivector[index[nv]][j] = atom->ivector[index[nv]][i];
-    else if (style[nv] == DVEC)
+      atom->ivector[index[nv]][i] = 0;
+    } else if (style[nv] == DVEC) {
       atom->dvector[index[nv]][j] = atom->dvector[index[nv]][i];
-    else if (style[nv] == IARRAY) {
+      atom->dvector[index[nv]][i] = 0.0;
+    } else if (style[nv] == IARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++)
-	atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k];
+      for (k = 0; k < ncol; k++) {
+	    atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k];
+        atom->iarray[index[nv]][i][k] = 0;
+      }
     } else if (style[nv] == DARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++)
-	atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k];
+      for (k = 0; k < ncol; k++) {
+	    atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k];
+        atom->darray[index[nv]][i][k] = 0.0;
+      }
     }
   }
 }
@@ -696,16 +702,24 @@ int FixPropertyAtom::pack_exchange(int i, double *buf)
     if (style[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
     else if (style[nv] == CHARGE) buf[m++] = atom->q[i];
     else if (style[nv] == RMASS) buf[m++] = atom->rmass[i];
-    else if (style[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
-    else if (style[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
-    else if (style[nv] == IARRAY) {
+    else if (style[nv] == IVEC) {
+      buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
+      atom->ivector[index[nv]][i] = 0;
+    } else if (style[nv] == DVEC) {
+      buf[m++] = atom->dvector[index[nv]][i];
+      atom->dvector[index[nv]][i] = 0;
+    } else if (style[nv] == IARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++)
-	buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
+      for (k = 0; k < ncol; k++) {
+	    buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
+        atom->iarray[index[nv]][i][k] = 0;
+      }
     } else if (style[nv] == DARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++)
-	buf[m++] = atom->darray[index[nv]][i][k];
+      for (k = 0; k < ncol; k++) {
+	    buf[m++] = atom->darray[index[nv]][i][k];
+        atom->darray[index[nv]][i][k] = 0.0;
+      }
     }
   }
 

src/fix_store_local.cpp

@@ -0,0 +1,200 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/ Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "fix_store_local.h"
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "memory.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define DELTA 1000
+
+/* ---------------------------------------------------------------------- */
+
+FixStoreLocal::FixStoreLocal(LAMMPS *lmp, int narg, char **arg) :
+    Fix(lmp, narg, arg), nvalues(0), vector(nullptr), array(nullptr)
+{
+  if (narg != 5) error->all(FLERR, "Illegal fix store/local command");
+  local_flag = 1;
+
+  nevery = utils::inumeric(FLERR, arg[3], false, lmp);
+  if (nevery <= 0) error->all(FLERR, "Illegal fix store/local command");
+  local_freq = nevery;
+
+  nvalues = utils::inumeric(FLERR, arg[4], false, lmp);
+  
+  if (nvalues <= 0) error->all(FLERR, "Illegal fix store/local command");
+  if (nvalues == 1)
+    size_local_cols = 0;
+  else
+    size_local_cols = nvalues;
+  size_local_rows = 0;
+  
+  vector = nullptr;
+  array = nullptr;  
+  nmax = 0;
+  ncount = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixStoreLocal::~FixStoreLocal()
+{
+  memory->destroy(vector);
+  memory->destroy(array);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixStoreLocal::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixStoreLocal::add_data(double *input_data, int i, int j)
+{
+  int *mask = atom->mask;
+  if (!(mask[i] & groupbit)) return;
+  if (!(mask[j] & groupbit)) return;
+  
+  if (ncount >= nmax) reallocate(ncount);
+  
+  // fill vector or array with local values
+  if (nvalues == 1) {
+    vector[ncount] = input_data[0];
+  } else {
+    for (int i = 0; i < nvalues; i++) array[ncount][i] = input_data[i];
+  }
+
+  ncount += 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixStoreLocal::post_force(int /*vflag*/)
+{
+  if (update->ntimestep % nevery == 0) {
+    size_local_rows = ncount;
+    ncount = 0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixStoreLocal::reallocate(int n)
+{
+  // grow vector or array
+  while (nmax <= n) nmax += DELTA;
+
+  if (nvalues == 1) {
+    memory->grow(vector, nmax, "fix_store_local:vector");
+    vector_local = vector;
+  } else {
+    memory->grow(array, nmax, nvalues, "fix_store_local:array");
+    array_local = array;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   write global array to restart file 
+------------------------------------------------------------------------- */
+
+void FixStoreLocal::write_restart(FILE *fp)
+{
+  // fill rbuf with size and vec/array values
+
+  rbuf[0] = nmax;
+  rbuf[1] = nvalues;
+  if (nvalues == 1) memcpy(&rbuf[2],vector,sizeof(double)*nmax);
+  else memcpy(&rbuf[2],&array[0][0],sizeof(double)*nmax*nvalues);
+
+  int n = nmax*nvalues + 2;
+  if (comm->me == 0) {
+    int size = n * sizeof(double);
+    fwrite(&size,sizeof(int),1,fp);
+    fwrite(rbuf,sizeof(double),n,fp);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   use global array from restart file to restart the Fix
+------------------------------------------------------------------------- */
+
+void FixStoreLocal::restart(char *buf)
+{
+  // first 2 values in buf are vec/array sizes
+
+  double *dbuf = (double *) buf;
+  int nrow_restart = dbuf[0];
+  int ncol_restart = dbuf[1];
+
+  // if size of vec/array has changed,
+  //   means the restart file is setting size of vec or array and doing init
+  //   because caller did not know size at time this fix was instantiated
+  // reallocate vector or array accordingly
+
+  if (nmax != nrow_restart || nvalues != ncol_restart) {
+    memory->destroy(vector);
+    memory->destroy(array);
+    memory->destroy(rbuf);
+    vector = nullptr;
+    array = nullptr;
+
+    nmax = nrow_restart;
+    nvalues = ncol_restart;
+    if (nvalues == 1) memory->create(vector,nmax,"fix/store/local:vector");
+    else memory->create(array,nmax,nvalues,"fix/store/local:array");
+    memory->create(rbuf,nmax*nvalues+2,"fix/store:rbuf");
+  }
+
+  int n = nmax*nvalues;
+  if (nvalues == 1) memcpy(vector,&dbuf[2],n*sizeof(double));
+  else memcpy(&array[0][0],&dbuf[2],n*sizeof(double));
+}
+
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
+------------------------------------------------------------------------- */
+
+int FixStoreLocal::maxsize_restart()
+{
+  return nvalues+1;
+}
+
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
+------------------------------------------------------------------------- */
+
+int FixStoreLocal::size_restart(int /*nlocal*/)
+{
+  return nvalues+1;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double FixStoreLocal::memory_usage()
+{
+  double bytes = (double) nmax * (double) nvalues * sizeof(double);
+  return bytes;
+}
+

src/fix_store_local.h

@@ -1,85 +1,75 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   https://www.lammps.org/ Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#ifdef FIX_CLASS
-// clang-format off
-FixStyle(pair/tracker,FixPairTracker);
-// clang-format on
-#else
-
-#ifndef LMP_FIX_PAIR_TRACKING_H
-#define LMP_FIX_PAIR_TRACKING_H
-
-#include "fix.h"
-
-namespace LAMMPS_NS {
-
-class FixPairTracker : public Fix {
- public:
-  FixPairTracker(class LAMMPS *, int, char **);
-  ~FixPairTracker();
-  int setmask();
-  void init();
-  void post_force(int);
-  double memory_usage();
-  void lost_contact(int, int, double, double, double, double);
-
- private:
-  int nvalues, nmax;
-  int index_i, index_j;
-  double tmin, rmin, rsum, time_initial, nstep_initial;
-
-  double *vector;
-  double **array;
-  int **type_filter;
-
-  int ncount;
-
-  void reallocate(int);
-
-  typedef void (FixPairTracker::*FnPtrPack)(int);
-  FnPtrPack *pack_choice;    // ptrs to pack functions
-
-  void pack_id1(int);
-  void pack_id2(int);
-
-  void pack_time_created(int);
-  void pack_time_broken(int);
-  void pack_time_total(int);
-
-  void pack_x(int);
-  void pack_y(int);
-  void pack_z(int);
-
-  void pack_rmin(int);
-  void pack_rave(int);
-};
-
-}    // namespace LAMMPS_NS
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory.  Check the input script syntax and compare to the
-documentation for the command.  You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-E: Invalid keyword in fix pair/tracker command
-
-Self-explanatory.
-
-*/
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/ Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(store/local,FixStoreLocal);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_STORE_LOCAL_H
+#define LMP_FIX_STORE_LOCAL_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixStoreLocal : public Fix {
+ public:
+  FixStoreLocal(class LAMMPS *, int, char **);
+  ~FixStoreLocal();
+  int setmask();
+  void post_force(int);
+  void write_restart(FILE *);
+  void restart(char *);
+  int size_restart(int);
+  int maxsize_restart();
+  double memory_usage();
+  void add_data(double *, int, int);
+  int nvalues;
+
+ private:
+  int nmax;
+
+  double *vector;
+  double **array;
+
+  int ncount;
+
+  void reallocate(int);
+  double *rbuf;    // restart buffer for GLOBAL vec/array
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Invalid keyword in fix store/local command
+
+Self-explanatory.
+
+E: Unused instance of fix store/local
+
+Instance of fix store/local is not associated with any other LAMMPS 
+class such as a bond style, pair style, etc.
+
+*/

src/lmptype.h

@@ -56,10 +56,12 @@
 namespace LAMMPS_NS {
 
 // reserve 2 hi bits in molecular system neigh list for special bonds flag
-// max local + ghost atoms per processor = 2^30 - 1
+// reserve 3rd last bit in neigh list for fix neigh/history flag
+// max local + ghost atoms per processor = 2^29 - 1
 
 #define SBBITS 30
-#define NEIGHMASK 0x3FFFFFFF
+#define HISTBITS 29
+#define NEIGHMASK 0x1FFFFFFF
 
 // default to 32-bit smallint and other ints, 64-bit bigint
 

src/math_extra.cpp

@@ -143,6 +143,58 @@ void richardson(double *q, double *m, double *w, double *moments, double dtq)
   MathExtra::qnormalize(q);
 }
 
+/* ----------------------------------------------------------------------
+   Richardson iteration to update quaternion from angular velocity
+   return new normalized quaternion q
+   also returns updated omega at 1/2 step
+   Assumes spherical particles - no need to rotate to match moments
+------------------------------------------------------------------------- */
+
+void richardson_sphere(double *q, double *w, double dtq)
+{
+  // full update from dq/dt = 1/2 w q
+
+  double wq[4];
+  MathExtra::vecquat(w,q,wq);
+
+  double qfull[4];
+  qfull[0] = q[0] + dtq * wq[0];
+  qfull[1] = q[1] + dtq * wq[1];
+  qfull[2] = q[2] + dtq * wq[2];
+  qfull[3] = q[3] + dtq * wq[3];
+  MathExtra::qnormalize(qfull);
+
+  // 1st half update from dq/dt = 1/2 w q
+
+  double qhalf[4];
+  qhalf[0] = q[0] + 0.5*dtq * wq[0];
+  qhalf[1] = q[1] + 0.5*dtq * wq[1];
+  qhalf[2] = q[2] + 0.5*dtq * wq[2];
+  qhalf[3] = q[3] + 0.5*dtq * wq[3];
+  MathExtra::qnormalize(qhalf);
+
+  // re-compute q at 1/2 step
+  // recompute wq
+
+  MathExtra::vecquat(w,qhalf,wq);
+
+  // 2nd half update from dq/dt = 1/2 w q
+
+  qhalf[0] += 0.5*dtq * wq[0];
+  qhalf[1] += 0.5*dtq * wq[1];
+  qhalf[2] += 0.5*dtq * wq[2];
+  qhalf[3] += 0.5*dtq * wq[3];
+  MathExtra::qnormalize(qhalf);
+
+  // corrected Richardson update
+
+  q[0] = 2.0*qhalf[0] - qfull[0];
+  q[1] = 2.0*qhalf[1] - qfull[1];
+  q[2] = 2.0*qhalf[2] - qfull[2];
+  q[3] = 2.0*qhalf[3] - qfull[3];
+  MathExtra::qnormalize(q);
+}
+
 /* ----------------------------------------------------------------------
    apply evolution operators to quat, quat momentum
    Miller et al., J Chem Phys. 116, 8649-8659 (2002)

src/math_extra.h

@@ -76,6 +76,7 @@ void write3(const double mat[3][3]);
 int mldivide3(const double mat[3][3], const double *vec, double *ans);
 void rotate(double matrix[3][3], int i, int j, int k, int l, double s, double tau);
 void richardson(double *q, double *m, double *w, double *moments, double dtq);
+void richardson_sphere(double *q, double *w, double dtq);
 void no_squish_rotate(int k, double *p, double *q, double *inertia, double dt);
 
 // shape matrix operations
@@ -91,6 +92,7 @@ inline void vecquat(double *a, double *b, double *c);
 inline void quatvec(double *a, double *b, double *c);
 inline void quatquat(double *a, double *b, double *c);
 inline void invquatvec(double *a, double *b, double *c);
+inline void quatrotvec(double *a, double *b, double *c);
 inline void axisangle_to_quat(const double *v, const double angle, double *quat);
 
 void angmom_to_omega(double *m, double *ex, double *ey, double *ez, double *idiag, double *w);
@@ -651,6 +653,29 @@ inline void MathExtra::invquatvec(double *a, double *b, double *c)
   c[2] = -a[3] * b[0] + a[2] * b[1] - a[1] * b[2] + a[0] * b[3];
 }
 
+/* ----------------------------------------------------------------------
+   quaternion rotation of vector: c = a*b*conj(a)
+   a is a quaternion
+   b is a three component vector
+   c is a three component vector
+------------------------------------------------------------------------- */
+
+inline void MathExtra::quatrotvec(double *a, double *b, double *c)
+{
+  double temp[4];
+
+  // temp = a*b
+  temp[0] = -a[1]*b[0] - a[2]*b[1] - a[3]*b[2];
+  temp[1] = a[0]*b[0] + a[2]*b[2] - a[3]*b[1];
+  temp[2] = a[0]*b[1] + a[3]*b[0] - a[1]*b[2];
+  temp[3] = a[0]*b[2] + a[1]*b[1] - a[2]*b[0];
+
+  // c = temp*conj(a)
+  c[0] = -a[1]*temp[0] + a[0]*temp[1] - a[3]*temp[2] + a[2]*temp[3];
+  c[1] = -a[2]*temp[0] + a[3]*temp[1] + a[0]*temp[2] - a[1]*temp[3];
+  c[2] = -a[3]*temp[0] - a[2]*temp[1] + a[1]*temp[2] + a[0]*temp[3];
+}
+
 /* ----------------------------------------------------------------------
    compute quaternion from axis-angle rotation
    v MUST be a unit vector

src/neighbor.cpp

@@ -21,6 +21,7 @@
 
 #include "atom.h"
 #include "atom_vec.h"
+#include "bond.h"
 #include "citeme.h"
 #include "comm.h"
 #include "compute.h"
@@ -1444,7 +1445,7 @@ void Neighbor::init_topology()
   // bonds,etc can only be broken for atom->molecular = Atom::MOLECULAR, not Atom::TEMPLATE
   // SHAKE sets bonds and angles negative
   // gcmc sets all bonds, angles, etc negative
-  // bond_quartic sets bonds to 0
+  // partial_flag sets bonds to 0
   // delete_bonds sets all interactions negative
 
   int bond_off = 0;
@@ -1453,7 +1454,9 @@ void Neighbor::init_topology()
     if (utils::strmatch(modify->fix[i]->style,"^shake")
         || utils::strmatch(modify->fix[i]->style,"^rattle"))
       bond_off = angle_off = 1;
-  if (force->bond && force->bond_match("quartic")) bond_off = 1;
+  if (force->bond)
+    if (force->bond->partial_flag)
+      bond_off = 1;
 
   if (atom->avec->bonds_allow && atom->molecular == Atom::MOLECULAR) {
     for (i = 0; i < atom->nlocal; i++) {

src/npair_half_size_bin_newtoff.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_bin_newtoff.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 
@@ -35,7 +38,8 @@ NPairHalfSizeBinNewtoff::NPairHalfSizeBinNewtoff(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeBinNewtoff::build(NeighList *list)
 {
-  int i,j,k,n,ibin;
+  int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -44,17 +48,26 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -68,6 +81,11 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
     ztmp = x[i][2];
     radi = radius[i];
     ibin = atom2bin[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in surrounding bins in stencil including self
     // only store pair if i < j
@@ -87,10 +105,24 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >= 0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/npair_half_size_bin_newton.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_bin_newton.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 
@@ -34,7 +37,8 @@ NPairHalfSizeBinNewton::NPairHalfSizeBinNewton(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeBinNewton::build(NeighList *list)
 {
-  int i,j,k,n,ibin;
+  int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -43,17 +47,26 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -66,6 +79,11 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over rest of atoms in i's bin, ghosts are at end of linked list
     // if j is owned atom, store it, since j is beyond i in linked list
@@ -90,10 +108,24 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >= 0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 
@@ -112,10 +144,24 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >= 0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/npair_half_size_bin_newton_tri.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_bin_newton_tri.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 
@@ -35,7 +38,8 @@ NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) :
 
 void NPairHalfSizeBinNewtonTri::build(NeighList *list)
 {
-  int i,j,k,n,ibin;
+  int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -44,17 +48,26 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -67,6 +80,11 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in bins in stencil
     // pairs for atoms j "below" i are excluded
@@ -96,10 +114,24 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >= 0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/npair_half_size_multi_newtoff.cpp

@@ -16,7 +16,10 @@ es   certain rights in this software.  This software is distributed under
 #include "npair_half_size_multi_newtoff.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -38,7 +41,9 @@ NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {
 
 void NPairHalfSizeMultiNewtoff::build(NeighList *list)
 {
-  int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
+  int which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -49,17 +54,26 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -73,6 +87,11 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     ibin = atom2bin[i];
 
@@ -93,27 +112,41 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list)
       ns = nstencil_multi[icollection][jcollection];
 
       for (k = 0; k < ns; k++) {
-            js = binhead_multi[jcollection][jbin + s[k]];
-            for (j = js; j >= 0; j = bins[j]) {
-              if (j <= i) continue;
+        js = binhead_multi[jcollection][jbin + s[k]];
+        for (j = js; j >= 0; j = bins[j]) {
+          if (j <= i) continue;
 
           jtype = type[j];
           if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
-              delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                  neighptr[n++] = j ^ mask_history;
-                else
-                  neighptr[n++] = j;
-              }
-            }
+          if (rsq <= cutdistsq) {
+            jh = j;
+            if (history && rsq < radsum*radsum)
+              jh = jh ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >= 0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = jh;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = jh;
+              else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+              // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+            } else neighptr[n++] = jh;
+          }
+        }
       }
     }
 

src/npair_half_size_multi_newton.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_multi_newton.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -36,7 +39,9 @@ NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeMultiNewton::build(NeighList *list)
 {
-  int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
+  int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
+  int which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -46,17 +51,26 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -70,6 +84,11 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     ibin = atom2bin[i];
 
@@ -94,33 +113,47 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
         //   if j is owned atom, store it if j > i
         //   if j is ghost, only store if j coords are "above and to the right" of i
 
-            for (j = js; j >= 0; j = bins[j]) {
+        for (j = js; j >= 0; j = bins[j]) {
           if ((icollection != jcollection) && (j < i)) continue;
 
-              if (j >= nlocal) {
-                if (x[j][2] < ztmp) continue;
-                if (x[j][2] == ztmp) {
-                  if (x[j][1] < ytmp) continue;
-                  if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
-                }
-              }
+          if (j >= nlocal) {
+            if (x[j][2] < ztmp) continue;
+            if (x[j][2] == ztmp) {
+              if (x[j][1] < ytmp) continue;
+              if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+            }
+          }
 
           jtype = type[j];
           if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
-              delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                  neighptr[n++] = j ^ mask_history;
-                else
-                  neighptr[n++] = j;
-              }
+          if (rsq <= cutdistsq) {
+            jh = j;
+            if (history && rsq < radsum*radsum)
+              jh = jh ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >= 0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = jh;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = jh;
+              else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+              // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+            } else neighptr[n++] = jh;
+          }
         }
       }
 
@@ -129,30 +162,44 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
       // stencil is half if i same size as j
       // stencil is full if i smaller than j
 
-          s = stencil_multi[icollection][jcollection];
-          ns = nstencil_multi[icollection][jcollection];
+      s = stencil_multi[icollection][jcollection];
+      ns = nstencil_multi[icollection][jcollection];
 
-          for (k = 0; k < ns; k++) {
-            js = binhead_multi[jcollection][jbin + s[k]];
-            for (j = js; j >= 0; j = bins[j]) {
+      for (k = 0; k < ns; k++) {
+        js = binhead_multi[jcollection][jbin + s[k]];
+        for (j = js; j >= 0; j = bins[j]) {
 
           jtype = type[j];
           if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
           delx = xtmp - x[j][0];
-              dely = ytmp - x[j][1];
-              delz = ztmp - x[j][2];
-              rsq = delx*delx + dely*dely + delz*delz;
-              radsum = radi + radius[j];
-              cutdistsq = (radsum+skin) * (radsum+skin);
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          radsum = radi + radius[j];
+          cutdistsq = (radsum+skin) * (radsum+skin);
 
-              if (rsq <= cutdistsq) {
-                if (history && rsq < radsum*radsum)
-                    neighptr[n++] = j ^ mask_history;
-                else
-                    neighptr[n++] = j;
-              }
-            }
+          if (rsq <= cutdistsq) {
+            jh = j;
+            if (history && rsq < radsum*radsum)
+              jh = jh ^ mask_history;
+
+            if (molecular != Atom::ATOMIC) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >= 0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = jh;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = jh;
+              else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+              // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+            } else neighptr[n++] = jh;
+          }
+        }
       }
     }
 

src/npair_half_size_multi_newton_tri.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_multi_newton_tri.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -36,7 +39,9 @@ NPairHalfSizeMultiNewtonTri::NPairHalfSizeMultiNewtonTri(LAMMPS *lmp) : NPair(lm
 
 void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
 {
-  int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
+  int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
+  int which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -46,17 +51,26 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -70,6 +84,11 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     ibin = atom2bin[i];
 
@@ -119,10 +138,24 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
               cutdistsq = (radsum+skin) * (radsum+skin);
 
               if (rsq <= cutdistsq) {
+                jh = j;
                 if (history && rsq < radsum*radsum)
-                    neighptr[n++] = j ^ mask_history;
-                else
-                    neighptr[n++] = j;
+                  jh = jh ^ mask_history;
+
+                if (molecular != Atom::ATOMIC) {
+                  if (!moltemplate)
+                    which = find_special(special[i],nspecial[i],tag[j]);
+                  else if (imol >= 0)
+                    which = find_special(onemols[imol]->special[iatom],
+                                         onemols[imol]->nspecial[iatom],
+                                         tag[j]-tagprev);
+                  else which = 0;
+                  if (which == 0) neighptr[n++] = jh;
+                  else if (domain->minimum_image_check(delx,dely,delz))
+                    neighptr[n++] = jh;
+                  else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+                  // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+                } else neighptr[n++] = jh;
               }
             }
           }

src/npair_half_size_multi_old_newtoff.cpp

@@ -15,8 +15,11 @@
 #include "npair_half_size_multi_old_newtoff.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
 #include "molecule.h"
+#include "molecule.h"
 #include "my_page.h"
 #include "neigh_list.h"
 
@@ -37,7 +40,8 @@ NPairHalfSizeMultiOldNewtoff::NPairHalfSizeMultiOldNewtoff(LAMMPS *lmp) : NPair(
 
 void NPairHalfSizeMultiOldNewtoff::build(NeighList *list)
 {
-  int i,j,k,n,itype,jtype,ibin,ns;
+  int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -47,17 +51,26 @@ void NPairHalfSizeMultiOldNewtoff::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -71,6 +84,11 @@ void NPairHalfSizeMultiOldNewtoff::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in other bins in stencil including self
     // only store pair if i < j
@@ -99,10 +117,24 @@ void NPairHalfSizeMultiOldNewtoff::build(NeighList *list)
         cutdistsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutdistsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >= 0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/npair_half_size_multi_old_newton.cpp

@@ -15,6 +15,8 @@
 #include "npair_half_size_multi_old_newton.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
 #include "molecule.h"
 #include "my_page.h"
@@ -36,7 +38,8 @@ NPairHalfSizeMultiOldNewton::NPairHalfSizeMultiOldNewton(LAMMPS *lmp) : NPair(lm
 
 void NPairHalfSizeMultiOldNewton::build(NeighList *list)
 {
-  int i,j,k,n,itype,jtype,ibin,ns;
+  int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -46,17 +49,26 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -70,6 +82,11 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over rest of atoms in i's bin, ghosts are at end of linked list
     // if j is owned atom, store it, since j is beyond i in linked list
@@ -95,10 +112,24 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list)
       cutdistsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutdistsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >= 0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 
@@ -125,10 +156,23 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list)
         cutdistsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutdistsq) {
-          if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
+        if (history && rsq < radsum*radsum)
+            j = j ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >= 0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = j;
+          else if (domain->minimum_image_check(delx,dely,delz))
             neighptr[n++] = j;
+          else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+          // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+        } else neighptr[n++] = j;
         }
       }
     }

src/npair_half_size_multi_old_newton_tri.cpp

@@ -15,6 +15,8 @@
 #include "npair_half_size_multi_old_newton_tri.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
 #include "molecule.h"
 #include "my_page.h"
@@ -35,7 +37,8 @@ NPairHalfSizeMultiOldNewtonTri::NPairHalfSizeMultiOldNewtonTri(LAMMPS *lmp) : NP
 
 void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
 {
-  int i,j,k,n,itype,jtype,ibin,ns;
+  int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutdistsq;
   int *neighptr,*s;
@@ -45,17 +48,26 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -69,7 +81,11 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
-
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over all atoms in bins, including self, in stencil
     // skip if i,j neighbor cutoff is less than bin distance
@@ -107,10 +123,24 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
         cutdistsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutdistsq) {
+          jh = j;
           if (history && rsq < radsum*radsum)
-            neighptr[n++] = j ^ mask_history;
-          else
-            neighptr[n++] = j;
+            jh = jh ^ mask_history;
+
+          if (molecular != Atom::ATOMIC) {
+            if (!moltemplate)
+              which = find_special(special[i],nspecial[i],tag[j]);
+            else if (imol >= 0)
+              which = find_special(onemols[imol]->special[iatom],
+                                   onemols[imol]->nspecial[iatom],
+                                   tag[j]-tagprev);
+            else which = 0;
+            if (which == 0) neighptr[n++] = jh;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = jh;
+            else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+            // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = jh;
         }
       }
     }

src/npair_half_size_nsq_newtoff.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_nsq_newtoff.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "group.h"
 #include "my_page.h"
 #include "neigh_list.h"
@@ -35,7 +38,8 @@ NPairHalfSizeNsqNewtoff::NPairHalfSizeNsqNewtoff(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeNsqNewtoff::build(NeighList *list)
 {
-  int i,j,n,bitmask;
+  int i,j,jh,n,bitmask,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -44,7 +48,10 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
+  tagint *tag = atom->tag;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   if (includegroup) {
@@ -52,13 +59,19 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
     bitmask = group->bitmask[includegroup];
   }
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -71,6 +84,11 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over remaining atoms, owned and ghost
 
@@ -86,10 +104,24 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >= 0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 

src/npair_half_size_nsq_newton.cpp

@@ -15,7 +15,10 @@
 #include "npair_half_size_nsq_newton.h"
 
 #include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
 #include "error.h"
+#include "molecule.h"
 #include "group.h"
 #include "my_page.h"
 #include "neigh_list.h"
@@ -36,7 +39,8 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeNsqNewton::build(NeighList *list)
 {
-  int i,j,n,itag,jtag,bitmask;
+  int i,j,jh,n,itag,jtag,bitmask,which,imol,iatom,moltemplate;
+  tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
   int *neighptr;
@@ -47,6 +51,8 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
   int *type = atom->type;
   int *mask = atom->mask;
   tagint *molecule = atom->molecule;
+  tagint **special = atom->special;
+  int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   if (includegroup) {
@@ -54,13 +60,19 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
     bitmask = group->bitmask[includegroup];
   }
 
+  int *molindex = atom->molindex;
+  int *molatom = atom->molatom;
+  Molecule **onemols = atom->avec->onemols;
+  if (molecular == Atom::TEMPLATE) moltemplate = 1;
+  else moltemplate = 0;
+
   int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  int mask_history = 3 << SBBITS;
+  int mask_history = 1 << HISTBITS;
 
   int inum = 0;
   ipage->reset();
@@ -74,6 +86,11 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
+    if (moltemplate) {
+      imol = molindex[i];
+      iatom = molatom[i];
+      tagprev = tag[i] - iatom - 1;
+    }
 
     // loop over remaining atoms, owned and ghost
 
@@ -105,10 +122,24 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
+        jh = j;
         if (history && rsq < radsum*radsum)
-          neighptr[n++] = j ^ mask_history;
-        else
-          neighptr[n++] = j;
+          jh = jh ^ mask_history;
+
+        if (molecular != Atom::ATOMIC) {
+          if (!moltemplate)
+            which = find_special(special[i],nspecial[i],tag[j]);
+          else if (imol >= 0)
+            which = find_special(onemols[imol]->special[iatom],
+                                 onemols[imol]->nspecial[iatom],
+                                 tag[j]-tagprev);
+          else which = 0;
+          if (which == 0) neighptr[n++] = jh;
+          else if (domain->minimum_image_check(delx,dely,delz))
+            neighptr[n++] = jh;
+          else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
+          // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
+        } else neighptr[n++] = jh;
       }
     }
 

src/set.cpp

@@ -303,7 +303,7 @@ void Set::command(int narg, char **arg)
       else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
       if (utils::strmatch(arg[iarg+4],"^v_")) varparse(arg[iarg+4],4);
       else wvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp);
-      if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag)
+      if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag)
         error->all(FLERR,"Cannot set this attribute for this atom style");
       set(QUAT);
       iarg += 5;
@@ -311,7 +311,7 @@ void Set::command(int narg, char **arg)
     } else if (strcmp(arg[iarg],"quat/random") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
-      if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag)
+      if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag)
         error->all(FLERR,"Cannot set this attribute for this atom style");
       if (ivalue <= 0)
         error->all(FLERR,"Invalid random number seed in set command");
@@ -937,18 +937,20 @@ void Set::set(int keyword)
       sp[i][3] = dvalue;
     }
 
-    // set quaternion orientation of ellipsoid or tri or body particle
-    // set quaternion orientation of ellipsoid or tri or body particle
+    // set quaternion orientation of ellipsoid or tri or body particle or sphere/bpm
     // enforce quat rotation vector in z dir for 2d systems
 
     else if (keyword == QUAT) {
       double *quat = nullptr;
+      double **quat2 = nullptr;
       if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
         quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
       else if (avec_tri && atom->tri[i] >= 0)
         quat = avec_tri->bonus[atom->tri[i]].quat;
       else if (avec_body && atom->body[i] >= 0)
         quat = avec_body->bonus[atom->body[i]].quat;
+      else if (atom->quat_flag)
+        quat2 = atom->quat;
       else
         error->one(FLERR,"Cannot set quaternion for atom that has none");
       if (domain->dimension == 2 && (xvalue != 0.0 || yvalue != 0.0))
@@ -957,11 +959,24 @@ void Set::set(int keyword)
 
       double theta2 = MY_PI2 * wvalue/180.0;
       double sintheta2 = sin(theta2);
-      quat[0] = cos(theta2);
-      quat[1] = xvalue * sintheta2;
-      quat[2] = yvalue * sintheta2;
-      quat[3] = zvalue * sintheta2;
-      MathExtra::qnormalize(quat);
+      if (atom->quat_flag) {
+        double temp[4];
+        temp[0] = cos(theta2);
+        temp[1] = xvalue * sintheta2;
+        temp[2] = yvalue * sintheta2;
+        temp[3] = zvalue * sintheta2;
+        MathExtra::qnormalize(temp);
+        quat2[i][0] = temp[0];
+        quat2[i][1] = temp[1];
+        quat2[i][2] = temp[2];
+        quat2[i][3] = temp[3];
+      } else {
+        quat[0] = cos(theta2);
+        quat[1] = xvalue * sintheta2;
+        quat[2] = yvalue * sintheta2;
+        quat[3] = zvalue * sintheta2;
+        MathExtra::qnormalize(quat);
+      }
     }
 
     // set theta of line particle
@@ -1223,6 +1238,7 @@ void Set::setrandom(int keyword)
   } else if (keyword == QUAT_RANDOM) {
     int nlocal = atom->nlocal;
     double *quat;
+    double **quat2;
 
     if (domain->dimension == 3) {
       double s,t1,t2,theta1,theta2;
@@ -1234,6 +1250,8 @@ void Set::setrandom(int keyword)
             quat = avec_tri->bonus[atom->tri[i]].quat;
           else if (avec_body && atom->body[i] >= 0)
             quat = avec_body->bonus[atom->body[i]].quat;
+          else if (atom->quat_flag)
+            quat2 = atom->quat;
           else
             error->one(FLERR,"Cannot set quaternion for atom that has none");
 
@@ -1243,10 +1261,17 @@ void Set::setrandom(int keyword)
           t2 = sqrt(s);
           theta1 = 2.0*MY_PI*ranpark->uniform();
           theta2 = 2.0*MY_PI*ranpark->uniform();
-          quat[0] = cos(theta2)*t2;
-          quat[1] = sin(theta1)*t1;
-          quat[2] = cos(theta1)*t1;
-          quat[3] = sin(theta2)*t2;
+          if (atom->quat_flag) {
+            quat2[i][0] = cos(theta2)*t2;
+            quat2[i][1] = sin(theta1)*t1;
+            quat2[i][2] = cos(theta1)*t1;
+            quat2[i][3] = sin(theta2)*t2;
+          } else {
+            quat[0] = cos(theta2)*t2;
+            quat[1] = sin(theta1)*t1;
+            quat[2] = cos(theta1)*t1;
+            quat[3] = sin(theta2)*t2;
+          }
           count++;
         }
 
@@ -1258,15 +1283,24 @@ void Set::setrandom(int keyword)
             quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
           else if (avec_body && atom->body[i] >= 0)
             quat = avec_body->bonus[atom->body[i]].quat;
+          else if (atom->quat_flag)
+            quat2 = atom->quat;
           else
             error->one(FLERR,"Cannot set quaternion for atom that has none");
 
           ranpark->reset(seed,x[i]);
           theta2 = MY_PI*ranpark->uniform();
-          quat[0] = cos(theta2);
-          quat[1] = 0.0;
-          quat[2] = 0.0;
-          quat[3] = sin(theta2);
+          if (atom->quat_flag) {
+            quat2[i][0] = cos(theta2);
+            quat2[i][1] = 0.0;
+            quat2[i][2] = 0.0;
+            quat2[i][3] = sin(theta2);
+          } else {
+            quat[0] = cos(theta2);
+            quat[1] = 0.0;
+            quat[2] = 0.0;
+            quat[3] = sin(theta2);
+          }
           count++;
         }
     }