diff --git a/src/BODY/Install.sh b/src/BODY/Install.sh
new file mode 100644
index 000000000..1836224e1
--- /dev/null
+++ b/src/BODY/Install.sh
@@ -0,0 +1,27 @@
+# Install/unInstall package files in LAMMPS
+
+if (test $1 = 1) then
+
+ cp body_nparticle.cpp ..
+ cp compute_body_local.cpp ..
+ cp fix_nve_body.cpp ..
+ cp pair_body.cpp ..
+
+ cp body_nparticle.h ..
+ cp compute_body_local.h ..
+ cp fix_nve_body.h ..
+ cp pair_body.h ..
+
+elif (test $1 = 0) then
+
+ rm -f ../body_nparticle.cpp
+ rm -f ../compute_body_local.cpp
+ rm -f ../fix_nve_body.cpp
+ rm -f ../pair_body.cpp
+
+ rm -f ../body_nparticle.h
+ rm -f ../compute_body_local.h
+ rm -f ../fix_nve_body.h
+ rm -f ../pair_body.h
+
+fi
diff --git a/src/BODY/body_nparticle.cpp b/src/BODY/body_nparticle.cpp
new file mode 100644
index 000000000..6d3ceecac
--- /dev/null
+++ b/src/BODY/body_nparticle.cpp
@@ -0,0 +1,201 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "body_nparticle.h"
+#include "math_extra.h"
+#include "atom_vec_body.h"
+#include "atom.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EPSILON 1.0e-7
+
+/* ---------------------------------------------------------------------- */
+
+BodyNparticle::BodyNparticle(LAMMPS *lmp, int narg, char **arg) :
+ Body(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Invalid body nparticle command");
+
+ int nmin = atoi(arg[1]);
+ int nmax = atoi(arg[2]);
+ if (nmin <= 0 || nmin > nmax)
+ error->all(FLERR,"Invalid body nparticle command");
+
+ size_forward = 0;
+ size_border = 1 + 3*nmax;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyNparticle::nsub(AtomVecBody::Bonus *bonus)
+{
+ return bonus->ivalue[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+double *BodyNparticle::coords(AtomVecBody::Bonus *bonus)
+{
+ return bonus->dvalue;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyNparticle::pack_border_body(AtomVecBody::Bonus *bonus, double *buf)
+{
+ int nsub = bonus->ivalue[0];
+ buf[0] = nsub;
+ memcpy(&buf[1],bonus->dvalue,3*nsub*sizeof(double));
+ return 1+3*nsub;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyNparticle::unpack_border_body(AtomVecBody::Bonus *bonus, double *buf)
+{
+ int nsub = static_cast<int> (buf[0]);
+ bonus->ivalue[0] = nsub;
+ memcpy(bonus->dvalue,&buf[1],3*nsub*sizeof(double));
+ return 1+3*nsub;
+}
+
+/* ----------------------------------------------------------------------
+ populate bonus data structure with data file values
+------------------------------------------------------------------------- */
+
+void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
+ char **ifile, char **dfile)
+{
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+
+ // error in data file if any values are NULL
+
+ for (int i = 0; i < ninteger; i++)
+ if (ifile[0] == NULL) error->one(FLERR,"");
+ for (int i = 0; i < ndouble; i++)
+ if (dfile[0] == NULL) error->one(FLERR,"");
+
+ // set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
+
+ if (ninteger != 1) error->one(FLERR,"");
+ int nsub = atoi(ifile[0]);
+ if (nsub < 1) error->one(FLERR,"");
+ if (ndouble != 6 + 3*nsub) error->one(FLERR,"");
+
+ bonus->ninteger = 1;
+ memory->create(bonus->ivalue,bonus->ninteger,"body:ivalue");
+ bonus->ivalue[0] = nsub;
+ bonus->ndouble = 3*nsub;
+ memory->create(bonus->dvalue,3*nsub,"body:dvalue");
+
+ // diagonalize inertia tensor
+
+ double tensor[3][3];
+ tensor[0][0] = atof(dfile[0]);
+ tensor[1][1] = atof(dfile[1]);
+ tensor[2][2] = atof(dfile[2]);
+ tensor[0][1] = tensor[1][0] = atof(dfile[3]);
+ tensor[0][2] = tensor[2][0] = atof(dfile[4]);
+ tensor[1][2] = tensor[2][1] = atof(dfile[5]);
+
+ double *inertia = bonus->inertia;
+ double evectors[3][3];
+ int ierror = MathExtra::jacobi(tensor,inertia,evectors);
+ if (ierror) error->one(FLERR,
+ "Insufficient Jacobi rotations for body nparticle");
+
+ // if any principal moment < scaled EPSILON, set to 0.0
+
+ double max;
+ max = MAX(inertia[0],inertia[1]);
+ max = MAX(max,inertia[2]);
+
+ if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
+ if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
+ if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
+
+ // exyz_space = principal axes in space frame
+
+ double ex_space[3],ey_space[3],ez_space[3];
+
+ ex_space[0] = evectors[0][0];
+ ex_space[1] = evectors[1][0];
+ ex_space[2] = evectors[2][0];
+ ey_space[0] = evectors[0][1];
+ ey_space[1] = evectors[1][1];
+ ey_space[2] = evectors[2][1];
+ ez_space[0] = evectors[0][2];
+ ez_space[1] = evectors[1][2];
+ ez_space[2] = evectors[2][2];
+
+ // enforce 3 evectors as a right-handed coordinate system
+ // flip 3rd vector if needed
+
+ double cross[3];
+ MathExtra::cross3(ex_space,ey_space,cross);
+ if (MathExtra::dot3(cross,ez_space) < 0.0) MathExtra::negate3(ez_space);
+
+ // create initial quaternion
+
+ MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus->quat);
+
+ // bonus->dvalue = sub-particle displacements in body frame
+
+ double delta[3],displace[3];
+
+ int j = 6;
+ int k = 0;
+ for (int i = 0; i < nsub; i++) {
+ delta[0] = atof(dfile[j]);
+ delta[1] = atof(dfile[j+1]);
+ delta[2] = atof(dfile[j+2]);
+ MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
+ delta,&bonus->dvalue[k]);
+ j += 3;
+ k += 3;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyNparticle::noutcol()
+{
+ return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyNparticle::noutrow(int ibonus)
+{
+ return avec->bonus[ibonus].ivalue[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BodyNparticle::output(int ibonus, int m, double *values)
+{
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*m],values);
+
+ double *x = atom->x[bonus->ilocal];
+ values[0] += x[0];
+ values[1] += x[1];
+ values[2] += x[2];
+}
diff --git a/src/BODY/body_nparticle.h b/src/BODY/body_nparticle.h
new file mode 100644
index 000000000..3263e3456
--- /dev/null
+++ b/src/BODY/body_nparticle.h
@@ -0,0 +1,57 @@
+/* ----------------------------------------------------------------------
+ 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 BODY_CLASS
+
+BodyStyle(nparticle,BodyNparticle)
+
+#else
+
+#ifndef LMP_BODY_NPARTICLE_H
+#define LMP_BODY_NPARTICLE_H
+
+#include "body.h"
+#include "atom_vec_body.h"
+
+namespace LAMMPS_NS {
+
+class BodyNparticle : public Body {
+ public:
+ BodyNparticle(class LAMMPS *, int, char **);
+ ~BodyNparticle() {}
+ int nsub(class AtomVecBody::Bonus *);
+ double *coords(class AtomVecBody::Bonus *);
+
+ int pack_border_body(class AtomVecBody::Bonus *, double *);
+ int unpack_border_body(class AtomVecBody::Bonus *, double *);
+ void data_body(int, int, int, char **, char **);
+
+ int noutrow(int);
+ int noutcol();
+ void output(int, int, double *);
+};
+
+}
+
+#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.
+
+*/
diff --git a/src/BODY/compute_body_local.cpp b/src/BODY/compute_body_local.cpp
new file mode 100644
index 000000000..a7a4fecfa
--- /dev/null
+++ b/src/BODY/compute_body_local.cpp
@@ -0,0 +1,206 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "compute_body_local.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "body.h"
+#include "update.h"
+#include "domain.h"
+#include "force.h"
+#include "bond.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define DELTA 10000
+
+enum{TYPE,INDEX};
+
+/* ---------------------------------------------------------------------- */
+
+ComputeBodyLocal::ComputeBodyLocal(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg < 4) error->all(FLERR,"Illegal compute body/local command");
+
+ local_flag = 1;
+ nvalues = narg - 3;
+ if (nvalues == 1) size_local_cols = 0;
+ else size_local_cols = nvalues;
+
+ which = new int[nvalues];
+ index = new int[nvalues];
+ nvalues = 0;
+
+ for (int iarg = 3; iarg < narg; iarg++) {
+ if (strcmp(arg[iarg],"type") == 0) which[nvalues++] = TYPE;
+ else {
+ which[nvalues] = INDEX;
+ index[nvalues] = force->inumeric(arg[iarg]) - 1;
+ nvalues++;
+ }
+ }
+
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec) error->all(FLERR,"Compute body/local requires atom style body");
+ bptr = avec->bptr;
+
+ int indexmax = bptr->noutcol();
+ for (int i = 0; i < nvalues; i++) {
+ if (which[i] == INDEX && (index[i] < 0 || index[i] >= indexmax))
+ error->all(FLERR,"Invalid index in compute body/local command");
+ }
+
+ nmax = 0;
+ vector = NULL;
+ array = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeBodyLocal::~ComputeBodyLocal()
+{
+ delete [] which;
+ delete [] index;
+ memory->destroy(vector);
+ memory->destroy(array);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeBodyLocal::init()
+{
+ // do initial memory allocation so that memory_usage() is correct
+
+ int ncount = compute_body(0);
+ if (ncount > nmax) reallocate(ncount);
+ size_local_rows = ncount;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeBodyLocal::compute_local()
+{
+ invoked_local = update->ntimestep;
+
+ // count local entries and compute body info
+
+ int ncount = compute_body(0);
+ if (ncount > nmax) reallocate(ncount);
+ size_local_rows = ncount;
+ ncount = compute_body(1);
+}
+
+/* ----------------------------------------------------------------------
+ count body info and compute body info on this proc
+ flag = 0, only count
+ flag = 1, also compute
+------------------------------------------------------------------------- */
+
+int ComputeBodyLocal::compute_body(int flag)
+{
+ // perform count
+
+ int *mask = atom->mask;
+ int *body = atom->body;
+ int nlocal = atom->nlocal;
+
+ int ncount = 0;
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ if (body[i] < 0) ncount++;
+ else ncount += bptr->noutrow(body[i]);
+ }
+
+ if (flag == 0) return ncount;
+
+ // perform computation and fill output vector/array
+
+ int m,n,ibonus;
+ double *values = new double[bptr->noutcol()];
+
+ double **x = atom->x;
+ int *type = atom->type;
+
+ ncount = 0;
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+ if (body[i] < 0) {
+ if (nvalues == 1) {
+ if (which[0] == TYPE) vector[ncount] = type[i];
+ else vector[ncount] = x[i][index[0]];
+ } else {
+ for (m = 0; m < nvalues; m++) {
+ if (which[m] == TYPE) array[ncount][m] = type[i];
+ else array[ncount][m] = x[i][index[m]];
+ }
+ }
+ ncount++;
+
+ } else {
+ ibonus = body[i];
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ n = bptr->noutrow(ibonus);
+ for (int j = 0; j < n; j++) {
+ bptr->output(ibonus,j,values);
+ if (nvalues == 1) {
+ if (which[0] == TYPE) vector[ncount] = type[i];
+ else vector[ncount] = values[index[0]];
+ } else {
+ for (m = 0; m < nvalues; m++) {
+ if (which[m] == TYPE) array[ncount][m] = type[i];
+ else array[ncount][m] = values[index[m]];
+ }
+ }
+ ncount++;
+ }
+ }
+ }
+ }
+
+ delete [] values;
+ return ncount;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeBodyLocal::reallocate(int n)
+{
+ // grow vector or array and indices array
+
+ while (nmax < n) nmax += DELTA;
+
+ if (nvalues == 1) {
+ memory->destroy(vector);
+ memory->create(vector,nmax,"body/local:vector");
+ vector_local = vector;
+ } else {
+ memory->destroy(array);
+ memory->create(array,nmax,nvalues,"body/local:array");
+ array_local = array;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeBodyLocal::memory_usage()
+{
+ double bytes = nmax*nvalues * sizeof(double);
+ return bytes;
+}
diff --git a/src/BODY/compute_body_local.h b/src/BODY/compute_body_local.h
new file mode 100644
index 000000000..6226c589f
--- /dev/null
+++ b/src/BODY/compute_body_local.h
@@ -0,0 +1,53 @@
+/* ----------------------------------------------------------------------
+ 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 COMPUTE_CLASS
+
+ComputeStyle(body/local,ComputeBodyLocal)
+
+#else
+
+#ifndef LMP_COMPUTE_BODY_LOCAL_H
+#define LMP_COMPUTE_BODY_LOCAL_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeBodyLocal : public Compute {
+ public:
+ ComputeBodyLocal(class LAMMPS *, int, char **);
+ ~ComputeBodyLocal();
+ void init();
+ void compute_local();
+ double memory_usage();
+
+ private:
+ int nvalues;
+ int *which,*index;
+
+ int nmax;
+ double *vector;
+ double **array;
+
+ class AtomVecBody *avec;
+ class Body *bptr;
+
+ int compute_body(int);
+ void reallocate(int);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/BODY/fix_nve_body.cpp b/src/BODY/fix_nve_body.cpp
new file mode 100644
index 000000000..27212ce0a
--- /dev/null
+++ b/src/BODY/fix_nve_body.cpp
@@ -0,0 +1,132 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "string.h"
+#include "fix_nve_body.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "force.h"
+#include "update.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEBody::FixNVEBody(LAMMPS *lmp, int narg, char **arg) :
+ FixNVE(lmp, narg, arg) {}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEBody::init()
+{
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec) error->all(FLERR,"Fix nve/body requires atom style body");
+
+ // check that all particles are bodies
+ // no point particles allowed
+
+ int *body = atom->body;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ if (body[i] < 0) error->one(FLERR,"Fix nve/body requires bodies");
+
+ FixNVE::init();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEBody::initial_integrate(int vflag)
+{
+ double dtfm;
+ double omega[3];
+ double *quat,*inertia;
+
+ AtomVecBody::Bonus *bonus = avec->bonus;
+ int *body = atom->body;
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double **angmom = atom->angmom;
+ double **torque = atom->torque;
+ 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;
+
+ 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];
+
+ // update angular momentum by 1/2 step
+
+ angmom[i][0] += dtf * torque[i][0];
+ angmom[i][1] += dtf * torque[i][1];
+ angmom[i][2] += dtf * torque[i][2];
+
+ // compute omega at 1/2 step from angmom at 1/2 step and current q
+ // update quaternion a full step via Richardson iteration
+ // returns new normalized quaternion
+
+ inertia = bonus[body[i]].inertia;
+ quat = bonus[body[i]].quat;
+ MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
+ MathExtra::richardson(quat,angmom[i],omega,inertia,dtq);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEBody::final_integrate()
+{
+ double dtfm;
+
+ double **v = atom->v;
+ double **f = atom->f;
+ double **angmom = atom->angmom;
+ double **torque = atom->torque;
+ double *rmass = atom->rmass;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+ 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];
+
+ angmom[i][0] += dtf * torque[i][0];
+ angmom[i][1] += dtf * torque[i][1];
+ angmom[i][2] += dtf * torque[i][2];
+ }
+}
diff --git a/src/BODY/fix_nve_body.h b/src/BODY/fix_nve_body.h
new file mode 100644
index 000000000..639e06d85
--- /dev/null
+++ b/src/BODY/fix_nve_body.h
@@ -0,0 +1,41 @@
+/* ----------------------------------------------------------------------
+ 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
+
+FixStyle(nve/body,FixNVEBody)
+
+#else
+
+#ifndef LMP_FIX_NVE_BODY_H
+#define LMP_FIX_NVE_BODY_H
+
+#include "fix_nve.h"
+
+namespace LAMMPS_NS {
+
+class FixNVEBody : public FixNVE {
+ public:
+ FixNVEBody(class LAMMPS *, int, char **);
+ void init();
+ void initial_integrate(int);
+ void final_integrate();
+
+ private:
+ double dtq;
+ class AtomVecBody *avec;
+};
+
+}
+#endif
+#endif
diff --git a/src/BODY/pair_body.cpp b/src/BODY/pair_body.cpp
new file mode 100644
index 000000000..5f9b4faea
--- /dev/null
+++ b/src/BODY/pair_body.cpp
@@ -0,0 +1,484 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_body.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "body_nparticle.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define DELTA 10000
+
+/* ---------------------------------------------------------------------- */
+
+PairBody::PairBody(LAMMPS *lmp) : Pair(lmp)
+{
+ dmax = nmax = 0;
+ discrete = NULL;
+ dnum = dfirst = NULL;
+
+ single_enable = 0;
+ restartinfo = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairBody::~PairBody()
+{
+ memory->destroy(discrete);
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(cut);
+ memory->destroy(epsilon);
+ memory->destroy(sigma);
+ memory->destroy(lj1);
+ memory->destroy(lj2);
+ memory->destroy(lj3);
+ memory->destroy(lj4);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBody::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ int ni,nj,npi,npj,ifirst,jfirst;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,r2inv,r6inv,forcelj;
+ double xi[3],xj[3],fi[3],fj[3],ti[3],tj[3];
+ double *dxi,*dxj;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double **torque = atom->torque;
+ int *body = atom->body;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int nall = nlocal + atom->nghost;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // grow discrete list if necessary and initialize
+
+ if (nall > nmax) {
+ nmax = nall;
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->create(dnum,nall,"pair:dnum");
+ memory->create(dfirst,nall,"pair:dfirst");
+ }
+ for (i = 0; i < nall; i++) dnum[i] = 0;
+ ndiscrete = 0;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq >= cutsq[itype][jtype]) continue;
+
+ // body/body interactions = NxM sub-particles
+
+ evdwl = 0.0;
+ if (body[i] >= 0 && body[j] >= 0) {
+ if (dnum[i] == 0) body2space(i);
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ if (dnum[j] == 0) body2space(j);
+ npj = dnum[j];
+ jfirst = dfirst[j];
+
+ for (ni = 0; ni < npi; ni++) {
+ dxi = discrete[ifirst+ni];
+
+ for (nj = 0; nj < npj; nj++) {
+ dxj = discrete[jfirst+nj];
+
+ xi[0] = x[i][0] + dxi[0];
+ xi[1] = x[i][1] + dxi[1];
+ xi[2] = x[i][2] + dxi[2];
+ xj[0] = x[j][0] + dxj[0];
+ xj[1] = x[j][1] + dxj[1];
+ xj[2] = x[j][2] + dxj[2];
+
+ delx = xi[0] - xj[0];
+ dely = xi[1] - xj[1];
+ delz = xi[2] - xj[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ fpair = forcelj*r2inv;
+
+ if (eflag)
+ evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+
+ fi[0] = delx*fpair;
+ fi[1] = dely*fpair;
+ fi[2] = delz*fpair;
+ f[i][0] += fi[0];
+ f[i][1] += fi[1];
+ f[i][2] += fi[2];
+ ti[0] = dxi[1]*fi[2] - dxi[2]*fi[1];
+ ti[1] = dxi[2]*fi[0] - dxi[0]*fi[2];
+ ti[2] = dxi[0]*fi[1] - dxi[1]*fi[0];
+ torque[i][0] += ti[0];
+ torque[i][1] += ti[1];
+ torque[i][2] += ti[2];
+
+ if (newton_pair || j < nlocal) {
+ fj[0] = -delx*fpair;
+ fj[1] = -dely*fpair;
+ fj[2] = -delz*fpair;
+ f[j][0] += fj[0];
+ f[j][1] += fj[1];
+ f[j][2] += fj[2];
+ tj[0] = dxj[1]*fj[2] - dxj[2]*fj[1];
+ tj[1] = dxj[2]*fj[0] - dxj[0]*fj[2];
+ tj[2] = dxj[0]*fj[1] - dxj[1]*fj[0];
+ torque[j][0] += tj[0];
+ torque[j][1] += tj[1];
+ torque[j][2] += tj[2];
+ }
+ }
+ }
+
+ // body/particle interaction = Nx1 sub-particles
+
+ } else if (body[i] >= 0) {
+ if (dnum[i] == 0) body2space(i);
+ npi = dnum[i];
+ ifirst = dfirst[i];
+
+ for (ni = 0; ni < npi; ni++) {
+ dxi = discrete[ifirst+ni];
+
+ xi[0] = x[i][0] + dxi[0];
+ xi[1] = x[i][1] + dxi[1];
+ xi[2] = x[i][2] + dxi[2];
+ xj[0] = x[j][0];
+ xj[1] = x[j][1];
+ xj[2] = x[j][2];
+
+ delx = xi[0] - xj[0];
+ dely = xi[1] - xj[1];
+ delz = xi[2] - xj[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ fpair = forcelj*r2inv;
+
+ if (eflag)
+ evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+
+ fi[0] = delx*fpair;
+ fi[1] = dely*fpair;
+ fi[2] = delz*fpair;
+ f[i][0] += fi[0];
+ f[i][1] += fi[1];
+ f[i][2] += fi[2];
+ ti[0] = dxi[1]*fi[2] - dxi[2]*fi[1];
+ ti[1] = dxi[2]*fi[0] - dxi[0]*fi[2];
+ ti[2] = dxi[0]*fi[1] - dxi[1]*fi[0];
+ torque[i][0] += ti[0];
+ torque[i][1] += ti[1];
+ torque[i][2] += ti[2];
+
+ if (newton_pair || j < nlocal) {
+ fj[0] = -delx*fpair;
+ fj[1] = -dely*fpair;
+ fj[2] = -delz*fpair;
+ f[j][0] += fj[0];
+ f[j][1] += fj[1];
+ f[j][2] += fj[2];
+ }
+ }
+
+
+ // particle/body interaction = Nx1 sub-particles
+
+ } else if (body[j] >= 0) {
+ if (dnum[j] == 0) body2space(j);
+ npj = dnum[j];
+ jfirst = dfirst[j];
+
+ for (nj = 0; nj < npj; nj++) {
+ dxj = discrete[jfirst+nj];
+
+ xi[0] = x[i][0];
+ xi[1] = x[i][1];
+ xi[2] = x[i][2];
+ xj[0] = x[j][0] + dxj[0];
+ xj[1] = x[j][1] + dxj[1];
+ xj[2] = x[j][2] + dxj[2];
+
+ delx = xi[0] - xj[0];
+ dely = xi[1] - xj[1];
+ delz = xi[2] - xj[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ fpair = forcelj*r2inv;
+
+ if (eflag)
+ evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+
+ fi[0] = delx*fpair;
+ fi[1] = dely*fpair;
+ fi[2] = delz*fpair;
+ f[i][0] += fi[0];
+ f[i][1] += fi[1];
+ f[i][2] += fi[2];
+
+ if (newton_pair || j < nlocal) {
+ fj[0] = -delx*fpair;
+ fj[1] = -dely*fpair;
+ fj[2] = -delz*fpair;
+ f[j][0] += fj[0];
+ f[j][1] += fj[1];
+ f[j][2] += fj[2];
+ tj[0] = dxj[1]*fj[2] - dxj[2]*fj[1];
+ tj[1] = dxj[2]*fj[0] - dxj[0]*fj[2];
+ tj[2] = dxj[0]*fj[1] - dxj[1]*fj[0];
+ torque[j][0] += tj[0];
+ torque[j][1] += tj[1];
+ torque[j][2] += tj[2];
+ }
+ }
+
+ // particle/particle interaction = 1x1 particles
+
+ } else {
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ fpair = forcelj*r2inv;
+
+ if (eflag)
+ evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairBody::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(cut,n+1,n+1,"pair:cut");
+ memory->create(epsilon,n+1,n+1,"pair:epsilon");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(lj1,n+1,n+1,"pair:lj1");
+ memory->create(lj2,n+1,n+1,"pair:lj2");
+ memory->create(lj3,n+1,n+1,"pair:lj3");
+ memory->create(lj4,n+1,n+1,"pair:lj4");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairBody::settings(int narg, char **arg)
+{
+ if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+
+ cut_global = force->numeric(arg[0]);
+
+ // reset cutoffs that have been explicitly set
+
+ if (allocated) {
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i+1; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairBody::coeff(int narg, char **arg)
+{
+ if (narg < 4 || narg > 5)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+ double epsilon_one = force->numeric(arg[2]);
+ double sigma_one = force->numeric(arg[3]);
+
+ double cut_one = cut_global;
+ if (narg == 5) cut_one = force->numeric(arg[4]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ epsilon[i][j] = epsilon_one;
+ sigma[i][j] = sigma_one;
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairBody::init_style()
+{
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec) error->all(FLERR,"Pair body requires atom style body");
+ if (strcmp(avec->bptr->style,"nparticle") != 0)
+ error->all(FLERR,"Pair body requires body style nparticle");
+ bptr = (BodyNparticle *) avec->bptr;
+
+ neighbor->request(this);
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairBody::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) {
+ epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+ sigma[i][i],sigma[j][j]);
+ sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+ cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+ }
+
+ lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+
+ epsilon[j][i] = epsilon[i][j];
+ sigma[j][i] = sigma[i][j];
+ lj1[j][i] = lj1[i][j];
+ lj2[j][i] = lj2[i][j];
+ lj3[j][i] = lj3[i][j];
+ lj4[j][i] = lj4[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ convert N sub-particles in body I to space frame using current quaternion
+ store sub-particle space-frame displacements from COM in discrete list
+------------------------------------------------------------------------- */
+
+void PairBody::body2space(int i)
+{
+ int ibonus = atom->body[i];
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int nsub = bptr->nsub(bonus);
+ double *coords = bptr->coords(bonus);
+
+ dnum[i] = nsub;
+ dfirst[i] = ndiscrete;
+
+ if (ndiscrete + nsub > dmax) {
+ dmax += DELTA;
+ memory->grow(discrete,dmax,3,"pair:discrete");
+ }
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+
+ for (int m = 0; m < nsub; m++) {
+ MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
+ ndiscrete++;
+ }
+}
diff --git a/src/BODY/pair_body.h b/src/BODY/pair_body.h
new file mode 100644
index 000000000..a1cc2cddf
--- /dev/null
+++ b/src/BODY/pair_body.h
@@ -0,0 +1,60 @@
+/* ----------------------------------------------------------------------
+ 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 PAIR_CLASS
+
+PairStyle(body,PairBody)
+
+#else
+
+#ifndef LMP_PAIR_BODY_H
+#define LMP_PAIR_BODY_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairBody : public Pair {
+ public:
+ PairBody(class LAMMPS *);
+ ~PairBody();
+ void compute(int, int);
+ void settings(int, char **);
+ void coeff(int, char **);
+ void init_style();
+ double init_one(int, int);
+
+ protected:
+ double cut_global;
+ double **cut;
+ double **epsilon,**sigma;
+ double **lj1,**lj2,**lj3,**lj4;
+
+ class AtomVecBody *avec;
+ class BodyNparticle *bptr;
+
+ double **discrete; // list of all sub-particles for all bodies
+ int ndiscrete; // number of discretes in list
+ int dmax; // allocated size of discrete list
+ int *dnum; // number of discretes per line, 0 if uninit
+ int *dfirst; // index of first discrete per each line
+ int nmax; // allocated size of dnum,dfirst vectors
+
+ void allocate();
+ void body2space(int);
+};
+
+}
+
+#endif
+#endif