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

src/DIPOLE/atom_vec_dipole.cpp

@@ -28,8 +28,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecDipole::AtomVecDipole(LAMMPS *lmp, int narg, char **arg) :
+AtomVecDipole::AtomVecDipole(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
   mass_type = 1;

src/DIPOLE/atom_vec_dipole.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecDipole : public AtomVec {
  public:
-  AtomVecDipole(class LAMMPS *, int, char **);
+  AtomVecDipole(class LAMMPS *);
   void grow(int);
   void grow_reset();
   void copy(int, int, int);

src/Depend.sh

@@ -46,4 +46,3 @@ elif (test $1 = 0) then
   fi
 
 fi
- 

src/MOLECULE/atom_vec_angle.cpp

@@ -27,8 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecAngle::AtomVecAngle(LAMMPS *lmp, int narg, char **arg) :
+AtomVecAngle::AtomVecAngle(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 1;
   bonds_allow = angles_allow = 1;

src/MOLECULE/atom_vec_angle.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecAngle : public AtomVec {
  public:
-  AtomVecAngle(class LAMMPS *, int, char **);
+  AtomVecAngle(class LAMMPS *);
   virtual ~AtomVecAngle() {}
   void grow(int);
   void grow_reset();

src/MOLECULE/atom_vec_bond.cpp

@@ -27,8 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecBond::AtomVecBond(LAMMPS *lmp, int narg, char **arg) :
+AtomVecBond::AtomVecBond(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 1;
   bonds_allow = 1;

src/MOLECULE/atom_vec_bond.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecBond : public AtomVec {
  public:
-  AtomVecBond(class LAMMPS *, int, char **);
+  AtomVecBond(class LAMMPS *);
   void grow(int);
   void grow_reset();
   void copy(int, int, int);

src/MOLECULE/atom_vec_full.cpp

@@ -27,8 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecFull::AtomVecFull(LAMMPS *lmp, int narg, char **arg) :
+AtomVecFull::AtomVecFull(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 1;
   bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 1;

src/MOLECULE/atom_vec_full.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecFull : public AtomVec {
  public:
-  AtomVecFull(class LAMMPS *, int, char **);
+  AtomVecFull(class LAMMPS *);
   virtual ~AtomVecFull() {}
   void grow(int);
   void grow_reset();

src/MOLECULE/atom_vec_molecular.cpp

@@ -27,8 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecMolecular::AtomVecMolecular(LAMMPS *lmp, int narg, char **arg) :
+AtomVecMolecular::AtomVecMolecular(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 1;
   bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 1;

src/MOLECULE/atom_vec_molecular.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecMolecular : public AtomVec {
  public:
-  AtomVecMolecular(class LAMMPS *, int, char **);
+  AtomVecMolecular(class LAMMPS *);
   void grow(int);
   void grow_reset();
   void copy(int, int, int);

src/Make.sh

@@ -10,7 +10,7 @@
 # else Make will not recreate them
 
 style () {
-  list=`grep -l $1 $2*.h`
+  list=`grep -sl $1 $2*.h`
   if (test -e style_$3.tmp) then
     rm -f style_$3.tmp
   fi
@@ -19,8 +19,15 @@ style () {
     echo "#include $qfile" >> style_$3.tmp
   done
   if (test ! -e style_$3.tmp) then
+    if (test "`cat style_$3.h`" != "") then
       rm -f style_$3.h
       touch style_$3.h
+      rm -f Obj_*/$4.d
+      if (test $5) then 
+        rm -f Obj_*/$5.d
+      fi
+      rm -f Obj_*/lammps.d
+      fi
   elif (test ! -e style_$3.h) then
     mv style_$3.tmp style_$3.h
     rm -f Obj_*/$4.d
@@ -50,7 +57,8 @@ style () {
 if (test $1 = "style") then
 
   style ANGLE_CLASS     angle_      angle      force
-  style ATOM_CLASS      atom_vec_   atom       atom
+  style ATOM_CLASS      atom_vec_   atom       atom      atom_vec_hybrid
+  style BODY_CLASS      body_       body       atom_vec_body
   style BOND_CLASS      bond_       bond       force
   style COMMAND_CLASS   ""          command    input
   style COMPUTE_CLASS   compute_    compute    modify

src/Makefile

@@ -13,7 +13,7 @@ OBJ = 	$(SRC:.cpp=.o)
 
 # Package variables
 
-PACKAGE = asphere class2 colloid dipole fld gpu granular kim \
+PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
 	  kspace manybody mc meam molecule opt peri poems reax replica \
 	  rigid shock srd voronoi xtc
 

src/PERI/atom_vec_peri.cpp

@@ -32,8 +32,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) :
+AtomVecPeri::AtomVecPeri(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
 

src/PERI/atom_vec_peri.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecPeri : public AtomVec {
  public:
-  AtomVecPeri(class LAMMPS *, int, char **);
+  AtomVecPeri(class LAMMPS *);
   void grow(int);
   void grow_reset();
   void copy(int, int, int);

src/RIGID/fix_rigid.cpp

@@ -2139,9 +2139,9 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
         array[id][0] = atof(values[5]);
         array[id][1] = atof(values[6]);
         array[id][2] = atof(values[7]);
-        array[id][5] = atof(values[8]);
-        array[id][4] = atof(values[9]);
         array[id][3] = atof(values[10]);
+        array[id][4] = atof(values[9]);
+        array[id][5] = atof(values[8]);
       }
 
       buf = next + 1;

src/USER-AWPMD/atom_vec_wavepacket.cpp

@@ -34,8 +34,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp, int narg, char **arg) :
+AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   comm_x_only = comm_f_only = 0;
 
@@ -56,8 +55,6 @@ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp, int narg, char **arg) :
     atom->ervel_flag = atom->erforce_flag = 1;
 
   atom->cs_flag = atom->csforce_flag = atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
-
-
 }
 
 /* ----------------------------------------------------------------------

src/USER-AWPMD/atom_vec_wavepacket.h

@@ -31,7 +31,7 @@ namespace LAMMPS_NS {
 
 class AtomVecWavepacket : public AtomVec {
 public:
-  AtomVecWavepacket(class LAMMPS *, int, char **);
+  AtomVecWavepacket(class LAMMPS *);
   ~AtomVecWavepacket() {}
   void grow(int);
   void grow_reset();

src/USER-CUDA/atom_vec_angle_cuda.cpp

@@ -60,8 +60,7 @@ using namespace LAMMPS_NS;
 #define BUF_FLOAT double
 /* ---------------------------------------------------------------------- */
 
-AtomVecAngleCuda::AtomVecAngleCuda(LAMMPS *lmp, int narg, char **arg) :
+AtomVecAngleCuda::AtomVecAngleCuda(LAMMPS *lmp) : AtomVecAngle(lmp)
-  AtomVecAngle(lmp, narg, arg)
 {
    cuda = lmp->cuda;
    if(cuda == NULL)

src/USER-CUDA/atom_vec_angle_cuda.h

@@ -37,7 +37,7 @@ namespace LAMMPS_NS {
 
 class AtomVecAngleCuda : public AtomVecAngle {
  public:
-  AtomVecAngleCuda(class LAMMPS *, int, char **);
+  AtomVecAngleCuda(class LAMMPS *);
   virtual ~AtomVecAngleCuda() {}
   void grow_copylist(int n);
   void grow_send(int n,double** buf_send,int flag);

src/USER-CUDA/atom_vec_atomic_cuda.cpp

@@ -60,8 +60,7 @@ using namespace LAMMPS_NS;
 #define BUF_FLOAT double
 /* ---------------------------------------------------------------------- */
 
-AtomVecAtomicCuda::AtomVecAtomicCuda(LAMMPS *lmp, int narg, char **arg) :
+AtomVecAtomicCuda::AtomVecAtomicCuda(LAMMPS *lmp) : AtomVecAtomic(lmp)
-  AtomVecAtomic(lmp, narg, arg)
 {
    cuda = lmp->cuda;
    if(cuda == NULL)

src/USER-CUDA/atom_vec_atomic_cuda.h

@@ -49,7 +49,7 @@ namespace LAMMPS_NS {
 
 class AtomVecAtomicCuda : public AtomVecAtomic {
  public:
-  AtomVecAtomicCuda(class LAMMPS *, int, char **);
+  AtomVecAtomicCuda(class LAMMPS *);
   virtual ~AtomVecAtomicCuda() {}
   void grow_copylist(int n);
   void grow_send(int n,double** buf_send,int flag);

src/USER-CUDA/atom_vec_charge_cuda.cpp

@@ -59,8 +59,7 @@ using namespace LAMMPS_NS;
 #define BUF_FLOAT double
 /* ---------------------------------------------------------------------- */
 
-AtomVecChargeCuda::AtomVecChargeCuda(LAMMPS *lmp, int narg, char **arg) :
+AtomVecChargeCuda::AtomVecChargeCuda(LAMMPS *lmp) : AtomVecCharge(lmp)
-  AtomVecCharge(lmp, narg, arg)
 {
    cuda = lmp->cuda;
    if(cuda == NULL)

src/USER-CUDA/atom_vec_charge_cuda.h

@@ -37,7 +37,7 @@ namespace LAMMPS_NS {
 
 class AtomVecChargeCuda : public AtomVecCharge {
  public:
-  AtomVecChargeCuda(class LAMMPS *, int, char **);
+  AtomVecChargeCuda(class LAMMPS *);
   virtual ~AtomVecChargeCuda() {}
   void grow_copylist(int n);
   void grow_send(int n,double** buf_send,int flag);

src/USER-EFF/atom_vec_electron.cpp

@@ -33,8 +33,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecElectron::AtomVecElectron(LAMMPS *lmp, int narg, char **arg) :
+AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   comm_x_only = comm_f_only = 0;
 

src/USER-EFF/atom_vec_electron.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecElectron : public AtomVec {
  public:
-  AtomVecElectron(class LAMMPS *, int, char **);
+  AtomVecElectron(class LAMMPS *);
   ~AtomVecElectron() {}
   void grow(int);
   void grow_reset();

src/USER-SPH/atom_vec_meso.cpp

@@ -27,8 +27,8 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecMeso::AtomVecMeso(LAMMPS *lmp, int narg, char **arg) :
+AtomVecMeso::AtomVecMeso(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg) {
+{
   molecular = 0;
   mass_type = 1;
 

src/USER-SPH/atom_vec_meso.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecMeso : public AtomVec {
  public:
-  AtomVecMeso(class LAMMPS *, int, char **);
+  AtomVecMeso(class LAMMPS *);
   ~AtomVecMeso() {}
   void grow(int);
   void grow_reset();

src/atom.cpp

@@ -42,6 +42,7 @@ using namespace LAMMPS_NS;
 #define DELTA_MEMSTR 1024
 #define EPSILON 1.0e-6
 #define CUDA_CHUNK 3000
+#define MAXBODY 20       // max # of lines in one body, also in ReadData class
 
 /* ---------------------------------------------------------------------- */
 
@@ -190,6 +191,7 @@ Atom::~Atom()
   memory->destroy(ellipsoid);
   memory->destroy(line);
   memory->destroy(tri);
+  memory->destroy(body);
   memory->destroy(spin);
   memory->destroy(eradius);
   memory->destroy(ervel);
@@ -271,7 +273,8 @@ void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
   vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
 
   int sflag;
-  avec = new_avec(style,narg,arg,suffix,sflag);
+  avec = new_avec(style,suffix,sflag);
+  avec->settings(narg,arg);
 
   if (sflag) {
     char estyle[256];
@@ -295,8 +298,7 @@ void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
    generate an AtomVec class, first with suffix appended
 ------------------------------------------------------------------------- */
 
-AtomVec *Atom::new_avec(const char *style, int narg, char **arg,
+AtomVec *Atom::new_avec(const char *style, char *suffix, int &sflag)
-                        char *suffix, int &sflag)
 {
   if (suffix && lmp->suffix_enable) {
     sflag = 1;
@@ -307,7 +309,7 @@ AtomVec *Atom::new_avec(const char *style, int narg, char **arg,
 
 #define ATOM_CLASS
 #define AtomStyle(key,Class) \
-    else if (strcmp(estyle,#key) == 0) return new Class(lmp,narg,arg);
+    else if (strcmp(estyle,#key) == 0) return new Class(lmp);
 #include "style_atom.h"
 #undef AtomStyle
 #undef ATOM_CLASS
@@ -320,7 +322,7 @@ AtomVec *Atom::new_avec(const char *style, int narg, char **arg,
 
 #define ATOM_CLASS
 #define AtomStyle(key,Class) \
-  else if (strcmp(style,#key) == 0) return new Class(lmp,narg,arg);
+  else if (strcmp(style,#key) == 0) return new Class(lmp);
 #include "style_atom.h"
 #undef ATOM_CLASS
 
@@ -718,6 +720,45 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus)
   delete [] values;
 }
 
+/* ----------------------------------------------------------------------
+   unpack n lines from atom-style specific section of data file
+   check that atom IDs are > 0 and <= map_tag_max
+   call style-specific routine to parse line
+------------------------------------------------------------------------- */
+
+void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body)
+{
+  int j,m,tagdata,ninteger,ndouble;
+
+  char **ivalues = new char*[10*MAXBODY];
+  char **dvalues = new char*[10*MAXBODY];
+
+  // loop over lines of body data
+  // tokenize the lines into ivalues and dvalues
+  // if I own atom tag, unpack its values
+
+  for (int i = 0; i < n; i++) {
+    if (i == 0) tagdata = atoi(strtok(buf," \t\n\r\f"));
+    else tagdata = atoi(strtok(NULL," \t\n\r\f"));
+    ninteger = atoi(strtok(NULL," \t\n\r\f"));
+    ndouble = atoi(strtok(NULL," \t\n\r\f"));
+
+    for (j = 0; j < ninteger; j++)
+      ivalues[j] = strtok(NULL," \t\n\r\f");
+    for (j = 0; j < ndouble; j++)
+      dvalues[j] = strtok(NULL," \t\n\r\f");
+
+    if (tagdata <= 0 || tagdata > map_tag_max)
+      error->one(FLERR,"Invalid atom ID in Bodies section of data file");
+
+    if ((m = map(tagdata)) >= 0)
+      avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);
+  }
+
+  delete [] ivalues;
+  delete [] dvalues;
+}
+
 /* ----------------------------------------------------------------------
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */

src/atom.h

@@ -53,7 +53,7 @@ class Atom : protected Pointers {
   double **omega,**angmom,**torque;
   double *radius,*rmass,*vfrac,*s0;
   double **x0;
-  int *ellipsoid,*line,*tri;
+  int *ellipsoid,*line,*tri,*body;
   int *spin;
   double *eradius,*ervel,*erforce,*ervelforce;
   double *cs,*csforce,*vforce;
@@ -89,7 +89,8 @@ class Atom : protected Pointers {
   // atom style and per-atom array existence flags
   // customize by adding new flag
 
-  int sphere_flag,ellipsoid_flag,line_flag,tri_flag,peri_flag,electron_flag;
+  int sphere_flag,ellipsoid_flag,line_flag,tri_flag,body_flag;
+  int peri_flag,electron_flag;
   int ecp_flag;
   int wavepacket_flag,sph_flag;
 
@@ -135,7 +136,7 @@ class Atom : protected Pointers {
 
   void settings(class Atom *);
   void create_avec(const char *, int, char **, char *suffix = NULL);
-  class AtomVec *new_avec(const char *, int, char **, char *, int &);
+  class AtomVec *new_avec(const char *, char *, int &);
   void init();
   void setup();
 
@@ -150,6 +151,7 @@ class Atom : protected Pointers {
   void data_atoms(int, char *);
   void data_vels(int, char *);
   void data_bonus(int, char *, class AtomVec *);
+  void data_bodies(int, char *, class AtomVecBody *);
 
   void data_bonds(int, char *);
   void data_angles(int, char *);

src/atom_vec.cpp

@@ -21,7 +21,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVec::AtomVec(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
+AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp)
 {
   nmax = 0;
   bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 0;
@@ -30,6 +30,15 @@ AtomVec::AtomVec(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   cudable = false;
 }
 
+/* ----------------------------------------------------------------------
+   no additional args by default
+------------------------------------------------------------------------- */
+
+void AtomVec::settings(int narg, char **arg)
+{
+  if (narg) error->all(FLERR,"Invalid atom_style command");
+}
+
 /* ----------------------------------------------------------------------
    copy of velocity remap settings from Domain
 ------------------------------------------------------------------------- */

src/atom_vec.h

@@ -14,6 +14,7 @@
 #ifndef LMP_ATOM_VEC_H
 #define LMP_ATOM_VEC_H
 
+#include "stdio.h"
 #include "pointers.h"
 
 namespace LAMMPS_NS {
@@ -41,8 +42,9 @@ class AtomVec : protected Pointers {
   int cudable;                         // 1 if atom style is CUDA-enabled
   int *maxsend;                        // CUDA-specific variable
 
-  AtomVec(class LAMMPS *, int, char **);
+  AtomVec(class LAMMPS *);
   virtual ~AtomVec() {}
+  virtual void settings(int, char **);
   virtual void init();
 
   virtual void grow(int) = 0;
@@ -76,6 +78,9 @@ class AtomVec : protected Pointers {
   virtual int pack_restart(int, double *) = 0;
   virtual int unpack_restart(double *) = 0;
 
+  virtual void write_restart_settings(FILE *) {}
+  virtual void read_restart_settings(FILE *) {}
+
   virtual void create_atom(int, double *) = 0;
   virtual void data_atom(double *, tagint, char **) = 0;
   virtual void data_atom_bonus(int, char **) {}

src/atom_vec_atomic.cpp

@@ -27,8 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecAtomic::AtomVecAtomic(LAMMPS *lmp, int narg, char **arg) :
+AtomVecAtomic::AtomVecAtomic(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
   mass_type = 1;

src/atom_vec_atomic.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecAtomic : public AtomVec {
  public:
-  AtomVecAtomic(class LAMMPS *, int, char **);
+  AtomVecAtomic(class LAMMPS *);
   virtual ~AtomVecAtomic() {}
   void grow(int);
   void grow_reset();

src/atom_vec_charge.cpp

@@ -27,8 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecCharge::AtomVecCharge(LAMMPS *lmp, int narg, char **arg) :
+AtomVecCharge::AtomVecCharge(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
   mass_type = 1;

src/atom_vec_charge.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecCharge : public AtomVec {
  public:
-  AtomVecCharge(class LAMMPS *, int, char **);
+  AtomVecCharge(class LAMMPS *);
   virtual ~AtomVecCharge() {}
   void grow(int);
   void grow_reset();

src/atom_vec_ellipsoid.cpp

@@ -36,8 +36,7 @@ using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) :
+AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
 
@@ -107,6 +106,7 @@ void AtomVecEllipsoid::grow_reset()
   mask = atom->mask; image = atom->image;
   x = atom->x; v = atom->v; f = atom->f;
   rmass = atom->rmass; angmom = atom->angmom; torque = atom->torque;
+  ellipsoid = atom->ellipsoid;
 }
 
 /* ----------------------------------------------------------------------
@@ -145,16 +145,16 @@ void AtomVecEllipsoid::copy(int i, int j, int delflag)
   angmom[j][1] = angmom[i][1];
   angmom[j][2] = angmom[i][2];
 
-  // if delflag and atom J has bonus data, then delete it
+  // if deleting atom J via delflag and J has bonus data, then delete it
 
   if (delflag && ellipsoid[j] >= 0) {
     copy_bonus(nlocal_bonus-1,ellipsoid[j]);
     nlocal_bonus--;
   }
 
-  // if atom I has bonus data and not deleting I, repoint I's bonus to J
+  // if atom I has bonus data, reset I's bonus.ilocal to loc J
 
-  if (ellipsoid[i] >= 0 && i != j) bonus[ellipsoid[i]].ilocal = j;
+  if (ellipsoid[i] >= 0) bonus[ellipsoid[i]].ilocal = j;
   ellipsoid[j] = ellipsoid[i];
 
   if (atom->nextra_grow)

src/atom_vec_ellipsoid.h

@@ -33,7 +33,7 @@ class AtomVecEllipsoid : public AtomVec {
   };
   struct Bonus *bonus;
 
-  AtomVecEllipsoid(class LAMMPS *, int, char **);
+  AtomVecEllipsoid(class LAMMPS *);
   ~AtomVecEllipsoid();
   void grow(int);
   void grow_reset();

src/atom_vec_hybrid.cpp

@@ -27,30 +27,62 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp, int narg, char **arg) :
+AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp) : AtomVec(lmp) {}
-  AtomVec(lmp, narg, arg)
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecHybrid::~AtomVecHybrid()
 {
-  int i,k,dummy;
+  for (int k = 0; k < nstyles; k++) delete styles[k];
-
+  delete [] styles;
-  if (narg < 1) error->all(FLERR,"Illegal atom_style command");
+  for (int k = 0; k < nstyles; k++) delete [] keywords[k];
-
+  delete [] keywords;
-  // create sub-styles
-
-  nstyles = narg;
-  styles = new AtomVec*[nstyles];
-  keywords = new char*[nstyles];
-
-  for (i = 0; i < narg; i++) {
-    for (k = 0; k < i; k++)
-      if (strcmp(arg[i],keywords[k]) == 0)
-        error->all(FLERR,"Atom style hybrid cannot use same atom style twice");
-    if (strcmp(arg[i],"hybrid") == 0)
-      error->all(FLERR,"Atom style hybrid cannot have hybrid as an argument");
-    styles[i] = atom->new_avec(arg[i],0,NULL,NULL,dummy);
-    keywords[i] = new char[strlen(arg[i])+1];
-    strcpy(keywords[i],arg[i]);
 }
 
+/* ----------------------------------------------------------------------
+   process sub-style args
+------------------------------------------------------------------------- */
+
+void AtomVecHybrid::settings(int narg, char **arg)
+{
+  // build list of all known atom styles
+
+  build_styles();
+
+  // allocate list of sub-styles as big as possibly needed if no extra args
+
+  styles = new AtomVec*[narg];
+  keywords = new char*[narg];
+
+  // allocate each sub-style
+  // call settings() with set of args that are not atom style names
+  // use known_style() to determine which args these are
+
+  int i,jarg,dummy;
+
+  int iarg = 0;
+  nstyles = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"hybrid") == 0)
+      error->all(FLERR,"Atom style hybrid cannot have hybrid as an argument");
+    for (i = 0; i < nstyles; i++)
+      if (strcmp(arg[iarg],keywords[i]) == 0)
+        error->all(FLERR,"Atom style hybrid cannot use same atom style twice");
+    styles[nstyles] = atom->new_avec(arg[iarg],NULL,dummy);
+    keywords[nstyles] = new char[strlen(arg[iarg])+1];
+    strcpy(keywords[nstyles],arg[iarg]);
+    jarg = iarg + 1;
+    while (jarg < narg && !known_style(arg[jarg])) jarg++;
+    styles[nstyles]->settings(jarg-iarg-1,&arg[iarg+1]);
+    iarg = jarg;
+    nstyles++;
+  }
+
+  // free allstyles created by build_styles()
+
+  for (int i = 0; i < nallstyles; i++) delete [] allstyles[i];
+  delete [] allstyles;
+
   // hybrid settings are MAX or MIN of sub-style settings
   // hybrid sizes are minimial values plus extra values for each sub-style
 
@@ -64,7 +96,7 @@ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp, int narg, char **arg) :
   size_data_vel = 4;
   xcol_data = 3;
 
-  for (k = 0; k < nstyles; k++) {
+  for (int k = 0; k < nstyles; k++) {
     molecular = MAX(molecular,styles[k]->molecular);
     bonds_allow = MAX(bonds_allow,styles[k]->bonds_allow);
     angles_allow = MAX(angles_allow,styles[k]->angles_allow);
@@ -89,16 +121,6 @@ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecHybrid::~AtomVecHybrid()
-{
-  for (int k = 0; k < nstyles; k++) delete styles[k];
-  delete [] styles;
-  for (int k = 0; k < nstyles; k++) delete [] keywords[k];
-  delete [] keywords;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void AtomVecHybrid::init()
 {
   AtomVec::init();
@@ -780,6 +802,22 @@ int AtomVecHybrid::unpack_restart(double *buf)
   return m;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void AtomVecHybrid::write_restart_settings(FILE *fp)
+{
+  for (int k = 0; k < nstyles; k++)
+    styles[k]->write_restart_settings(fp);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecHybrid::read_restart_settings(FILE *fp)
+{
+  for (int k = 0; k < nstyles; k++)
+    styles[k]->read_restart_settings(fp);
+}
+
 /* ----------------------------------------------------------------------
    create one atom of itype at coord
    create each sub-style one after the other
@@ -863,6 +901,45 @@ void AtomVecHybrid::data_vel(int m, char **values)
     n += styles[k]->data_vel_hybrid(m,&values[n]);
 }
 
+/* ----------------------------------------------------------------------
+   allstyles = list of all atom styles in this LAMMPS executable
+------------------------------------------------------------------------- */
+
+void AtomVecHybrid::build_styles()
+{
+  nallstyles = 0;
+#define ATOM_CLASS
+#define AtomStyle(key,Class) nallstyles++;
+#include "style_atom.h"
+#undef AtomStyle
+#undef ATOM_CLASS
+
+  allstyles = new char*[nallstyles];
+
+  int n;
+  nallstyles = 0;
+#define ATOM_CLASS
+#define AtomStyle(key,Class)                \
+  n = strlen(#key) + 1;                     \
+  allstyles[nallstyles] = new char[n];      \
+  strcpy(allstyles[nallstyles],#key);       \
+  nallstyles++;
+#include "style_atom.h"
+#undef AtomStyle
+#undef ATOM_CLASS
+}
+
+/* ----------------------------------------------------------------------
+   allstyles = list of all known atom styles
+------------------------------------------------------------------------- */
+
+int AtomVecHybrid::known_style(char *str)
+{
+  for (int i = 0; i < nallstyles; i++)
+    if (strcmp(str,allstyles[i]) == 0) return 1;
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
 ------------------------------------------------------------------------- */

src/atom_vec_hybrid.h

@@ -20,6 +20,7 @@ AtomStyle(hybrid,AtomVecHybrid)
 #ifndef LMP_ATOM_VEC_HYBRID_H
 #define LMP_ATOM_VEC_HYBRID_H
 
+#include "stdio.h"
 #include "atom_vec.h"
 
 namespace LAMMPS_NS {
@@ -30,8 +31,9 @@ class AtomVecHybrid : public AtomVec {
   class AtomVec **styles;
   char **keywords;
 
-  AtomVecHybrid(class LAMMPS *, int, char **);
+  AtomVecHybrid(class LAMMPS *);
   ~AtomVecHybrid();
+  void settings(int, char **);
   void init();
   void grow(int);
   void grow_reset();
@@ -52,6 +54,8 @@ class AtomVecHybrid : public AtomVec {
   int size_restart();
   int pack_restart(int, double *);
   int unpack_restart(double *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
   void create_atom(int, double *);
   void data_atom(double *, tagint, char **);
   int data_atom_hybrid(int, char **) {return 0;}
@@ -63,6 +67,12 @@ class AtomVecHybrid : public AtomVec {
   tagint *image;
   double **x,**v,**f;
   double **omega,**angmom;
+
+  int nallstyles;
+  char **allstyles;
+
+  void build_styles();
+  int known_style(char *);
 };
 
 }

src/atom_vec_line.cpp

@@ -32,8 +32,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecLine::AtomVecLine(LAMMPS *lmp, int narg, char **arg) :
+AtomVecLine::AtomVecLine(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
 
@@ -116,6 +115,7 @@ void AtomVecLine::grow_reset()
   x = atom->x; v = atom->v; f = atom->f;
   molecule = atom->molecule; rmass = atom->rmass;
   omega = atom->omega; torque = atom->torque;
+  line = atom->line;
 }
 
 /* ----------------------------------------------------------------------
@@ -155,16 +155,16 @@ void AtomVecLine::copy(int i, int j, int delflag)
   omega[j][1] = omega[i][1];
   omega[j][2] = omega[i][2];
 
-  // if delflag and atom J has bonus data, then delete it
+  // if deleting atom J via delflag and J has bonus data, then delete it
 
   if (delflag && line[j] >= 0) {
     copy_bonus(nlocal_bonus-1,line[j]);
     nlocal_bonus--;
   }
 
-  // if atom I has bonus data and not deleting I, repoint I's bonus to J
+  // if atom I has bonus data, reset I's bonus.ilocal to loc J
 
-  if (line[i] >= 0 && i != j) bonus[line[i]].ilocal = j;
+  if (line[i] >= 0) bonus[line[i]].ilocal = j;
   line[j] = line[i];
 
   if (atom->nextra_grow)

src/atom_vec_line.h

@@ -32,7 +32,7 @@ class AtomVecLine : public AtomVec {
   };
   struct Bonus *bonus;
 
-  AtomVecLine(class LAMMPS *, int, char **);
+  AtomVecLine(class LAMMPS *);
   ~AtomVecLine();
   void init();
   void grow(int);

src/atom_vec_sphere.cpp

@@ -33,8 +33,7 @@ using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecSphere::AtomVecSphere(LAMMPS *lmp, int narg, char **arg) :
+AtomVecSphere::AtomVecSphere(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
 

src/atom_vec_sphere.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class AtomVecSphere : public AtomVec {
  public:
-  AtomVecSphere(class LAMMPS *, int, char **);
+  AtomVecSphere(class LAMMPS *);
   ~AtomVecSphere() {}
   void init();
   void grow(int);

src/atom_vec_tri.cpp

@@ -33,8 +33,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecTri::AtomVecTri(LAMMPS *lmp, int narg, char **arg) :
+AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
-  AtomVec(lmp, narg, arg)
 {
   molecular = 0;
 
@@ -117,6 +116,7 @@ void AtomVecTri::grow_reset()
   x = atom->x; v = atom->v; f = atom->f;
   molecule = atom->molecule; rmass = atom->rmass;
   angmom = atom->angmom; torque = atom->torque;
+  tri = atom->tri;
 }
 
 /* ----------------------------------------------------------------------
@@ -157,16 +157,16 @@ void AtomVecTri::copy(int i, int j, int delflag)
   angmom[j][1] = angmom[i][1];
   angmom[j][2] = angmom[i][2];
 
-  // if delflag and atom J has bonus data, then delete it
+  // if deleting atom J via delflag and J has bonus data, then delete it
 
   if (delflag && tri[j] >= 0) {
     copy_bonus(nlocal_bonus-1,tri[j]);
     nlocal_bonus--;
   }
 
-  // if atom I has bonus data and not deleting I, repoint I's bonus to J
+  // if atom I has bonus data, reset I's bonus.ilocal to loc J
 
-  if (tri[i] >= 0 && i != j) bonus[tri[i]].ilocal = j;
+  if (tri[i] >= 0) bonus[tri[i]].ilocal = j;
   tri[j] = tri[i];
 
   if (atom->nextra_grow)
@@ -1320,7 +1320,7 @@ void AtomVecTri::create_atom(int itype, double *coord)
 }
 
 /* ----------------------------------------------------------------------
-   unpack one tri from Atoms section of data file
+   unpack one line from Atoms section of data file
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
@@ -1387,7 +1387,7 @@ int AtomVecTri::data_atom_hybrid(int nlocal, char **values)
 }
 
 /* ----------------------------------------------------------------------
-   unpack one tri from Tris section of data file
+   unpack one line from Tris section of data file
 ------------------------------------------------------------------------- */
 
 void AtomVecTri::data_atom_bonus(int m, char **values)
@@ -1511,7 +1511,7 @@ void AtomVecTri::data_atom_bonus(int m, char **values)
 }
 
 /* ----------------------------------------------------------------------
-   unpack one tri from Velocities section of data file
+   unpack one line from Velocities section of data file
 ------------------------------------------------------------------------- */
 
 void AtomVecTri::data_vel(int m, char **values)
@@ -1555,7 +1555,8 @@ bigint AtomVecTri::memory_usage()
   if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
   if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
   if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
-  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax*comm->nthreads,3);
+  if (atom->memcheck("torque")) bytes += 
+                                  memory->usage(torque,nmax*comm->nthreads,3);
   if (atom->memcheck("tri")) bytes += memory->usage(tri,nmax);
 
   bytes += nmax_bonus*sizeof(Bonus);

src/atom_vec_tri.h

@@ -34,7 +34,7 @@ class AtomVecTri : public AtomVec {
   };
   struct Bonus *bonus;
 
-  AtomVecTri(class LAMMPS *, int, char **);
+  AtomVecTri(class LAMMPS *);
   ~AtomVecTri();
   void init();
   void grow(int);

src/comm.cpp

@@ -800,8 +800,11 @@ void Comm::exchange()
   // b/c atoms migrate to new procs in exchange() and
   //   new ghosts are created in borders()
   // map_set() is done at end of borders()
+  // clear ghost count and any ghost bonus data internal to AtomVec
 
   if (map_style) atom->map_clear();
+  atom->nghost = 0;
+  atom->avec->clear_bonus();
 
   // subbox bounds for orthogonal or triclinic
 
@@ -836,6 +839,7 @@ void Comm::exchange()
     }
     atom->nlocal = nlocal;
 
+
     // send/recv atoms in both directions
     // if 1 proc in dimension, no send/recv, set recv buf to send buf
     // if 2 procs in dimension, single send/recv
@@ -907,11 +911,6 @@ void Comm::borders()
   MPI_Status status;
   AtomVec *avec = atom->avec;
 
-  // clear old ghosts and any ghost bonus data internal to AtomVec
-
-  atom->nghost = 0;
-  atom->avec->clear_bonus();
-
   // do swaps over all 3 dimensions
 
   iswap = 0;

src/compute_property_atom.cpp

@@ -19,6 +19,7 @@
 #include "atom_vec_ellipsoid.h"
 #include "atom_vec_line.h"
 #include "atom_vec_tri.h"
+#include "atom_vec_body.h"
 #include "update.h"
 #include "domain.h"
 #include "memory.h"
@@ -191,26 +192,36 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
       if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
                                       "atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_shapez;
+
     } else if (strcmp(arg[iarg],"quatw") == 0) {
       avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
-      if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+      avec_body = (AtomVecBody *) atom->style_match("body");
+      if (!avec_ellipsoid && !avec_body) 
+        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");
-      if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+      avec_body = (AtomVecBody *) atom->style_match("body");
+      if (!avec_ellipsoid && !avec_body) 
+        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");
-      if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+      avec_body = (AtomVecBody *) atom->style_match("body");
+      if (!avec_ellipsoid && !avec_body) 
+        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");
-      if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+      avec_body = (AtomVecBody *) atom->style_match("body");
+      if (!avec_ellipsoid && !avec_body) 
+        error->all(FLERR,"Compute property/atom for "
                    "atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_quatk;
+
     } else if (strcmp(arg[iarg],"tqx") == 0) {
       if (!atom->torque_flag)
         error->all(FLERR,"Compute property/atom for "
@@ -349,6 +360,7 @@ void ComputePropertyAtom::init()
   avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
   avec_line = (AtomVecLine *) atom->style_match("line");
   avec_tri = (AtomVecTri *) atom->style_match("tri");
+  avec_body = (AtomVecBody *) atom->style_match("body");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1136,6 +1148,7 @@ void ComputePropertyAtom::pack_shapez(int n)
 
 void ComputePropertyAtom::pack_quatw(int n)
 {
+  if (avec_ellipsoid) {
     AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
     int *ellipsoid = atom->ellipsoid;
     int *mask = atom->mask;
@@ -1147,12 +1160,27 @@ void ComputePropertyAtom::pack_quatw(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+
+  } else {
+    AtomVecBody::Bonus *bonus = avec_body->bonus;
+    int *body = atom->body;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if ((mask[i] & groupbit) && body[i] >= 0)
+        buf[n] = bonus[body[i]].quat[0];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyAtom::pack_quati(int n)
 {
+  if (avec_ellipsoid) {
     AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
     int *ellipsoid = atom->ellipsoid;
     int *mask = atom->mask;
@@ -1164,12 +1192,27 @@ void ComputePropertyAtom::pack_quati(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+
+  } else {
+    AtomVecBody::Bonus *bonus = avec_body->bonus;
+    int *body = atom->body;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if ((mask[i] & groupbit) && body[i] >= 0)
+        buf[n] = bonus[body[i]].quat[1];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyAtom::pack_quatj(int n)
 {
+  if (avec_ellipsoid) {
     AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
     int *ellipsoid = atom->ellipsoid;
     int *mask = atom->mask;
@@ -1181,12 +1224,27 @@ void ComputePropertyAtom::pack_quatj(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+
+  } else {
+    AtomVecBody::Bonus *bonus = avec_body->bonus;
+    int *body = atom->body;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if ((mask[i] & groupbit) && body[i] >= 0)
+        buf[n] = bonus[body[i]].quat[2];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyAtom::pack_quatk(int n)
 {
+  if (avec_ellipsoid) {
     AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
     int *ellipsoid = atom->ellipsoid;
     int *mask = atom->mask;
@@ -1198,6 +1256,20 @@ void ComputePropertyAtom::pack_quatk(int n)
       else buf[n] = 0.0;
       n += nvalues;
     }
+
+  } else {
+    AtomVecBody::Bonus *bonus = avec_body->bonus;
+    int *body = atom->body;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++) {
+      if ((mask[i] & groupbit) && body[i] >= 0)
+        buf[n] = bonus[body[i]].quat[3];
+      else buf[n] = 0.0;
+      n += nvalues;
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/compute_property_atom.h

@@ -41,6 +41,7 @@ class ComputePropertyAtom : public Compute {
   class AtomVecEllipsoid *avec_ellipsoid;
   class AtomVecLine *avec_line;
   class AtomVecTri *avec_tri;
+  class AtomVecBody *avec_body;
 
   typedef void (ComputePropertyAtom::*FnPtrPack)(int);
   FnPtrPack *pack_choice;              // ptrs to pack functions

src/dump_local.cpp

@@ -194,6 +194,17 @@ void DumpLocal::write_header(bigint ndump)
     fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep);
     fprintf(fp,"ITEM: NUMBER OF %s\n",label);
     fprintf(fp,BIGINT_FORMAT "\n",ndump);
+    if (domain->triclinic) {
+      fprintf(fp,"ITEM: BOX BOUNDS xy xz yz %s\n",boundstr);
+      fprintf(fp,"%g %g %g\n",boxxlo,boxxhi,boxxy);
+      fprintf(fp,"%g %g %g\n",boxylo,boxyhi,boxxz);
+      fprintf(fp,"%g %g %g\n",boxzlo,boxzhi,boxyz);
+    } else {
+      fprintf(fp,"ITEM: BOX BOUNDS %s\n",boundstr);
+      fprintf(fp,"%g %g\n",boxxlo,boxxhi);
+      fprintf(fp,"%g %g\n",boxylo,boxyhi);
+      fprintf(fp,"%g %g\n",boxzlo,boxzhi);
+    }
     fprintf(fp,"ITEM: %s %s\n",label,columns);
   }
 }
@@ -223,13 +234,15 @@ int DumpLocal::count()
   for (int i = 0; i < ncompute; i++) {
     if (nmine < 0) nmine = compute[i]->size_local_rows;
     else if (nmine != compute[i]->size_local_rows)
-      error->one(FLERR,"Dump local count is not consistent across input fields");
+      error->one(FLERR,
+                 "Dump local count is not consistent across input fields");
   }
 
   for (int i = 0; i < nfix; i++) {
     if (nmine < 0) nmine = fix[i]->size_local_rows;
     else if (nmine != fix[i]->size_local_rows)
-      error->one(FLERR,"Dump local count is not consistent across input fields");
+      error->one(FLERR,
+                 "Dump local count is not consistent across input fields");
   }
 
   return nmine;

src/read_data.cpp

@@ -44,9 +44,10 @@ using namespace LAMMPS_NS;
 #define LB_FACTOR 1.1
 #define CHUNK 1024
 #define DELTA 4            // must be 2 or larger
+#define MAXBODY 20         // max # of lines in one body, also in Atom class
 
                            // customize for new sections
-#define NSECTIONS 23       // change when add to header::section_keywords
+#define NSECTIONS 24       // change when add to header::section_keywords
 
 /* ---------------------------------------------------------------------- */
 
@@ -68,6 +69,8 @@ ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp)
   avec_line = (AtomVecLine *) atom->style_match("line");
   ntris = 0;
   avec_tri = (AtomVecTri *) atom->style_match("tri");
+  nbodies = 0;
+  avec_body = (AtomVecBody *) atom->style_match("body");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -219,6 +222,11 @@ void ReadData::command(int narg, char **arg)
         error->all(FLERR,"Invalid data file section: Triangles");
       if (atomflag == 0) error->all(FLERR,"Must read Atoms before Triangles");
       bonus(ntris,(AtomVec *) avec_tri,"triangles");
+    } else if (strcmp(keyword,"Bodies") == 0) {
+      if (!avec_body)
+        error->all(FLERR,"Invalid data file section: Bodies");
+      if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bodies");
+      bodies();
 
     } else if (strcmp(keyword,"Bonds") == 0) {
       if (atom->avec->bonds_allow == 0)
@@ -377,7 +385,7 @@ void ReadData::header(int flag)
   // customize for new sections
 
   const char *section_keywords[NSECTIONS] =
-    {"Atoms","Velocities","Ellipsoids","Lines","Triangles",
+    {"Atoms","Velocities","Ellipsoids","Lines","Triangles","Bodies",
      "Bonds","Angles","Dihedrals","Impropers",
      "Masses","Pair Coeffs","Bond Coeffs","Angle Coeffs",
      "Dihedral Coeffs","Improper Coeffs",
@@ -452,6 +460,10 @@ void ReadData::header(int flag)
       if (!avec_tri && me == 0)
         error->one(FLERR,"No triangles allowed with this atom style");
       sscanf(line,BIGINT_FORMAT,&ntris);
+    } else if (strstr(line,"bodies")) {
+      if (!avec_body && me == 0)
+        error->one(FLERR,"No bodies allowed with this atom style");
+      sscanf(line,BIGINT_FORMAT,&nbodies);
     }
 
     else if (strstr(line,"bonds")) sscanf(line,BIGINT_FORMAT,&atom->nbonds);
@@ -718,6 +730,86 @@ void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type)
   }
 }
 
+/* ----------------------------------------------------------------------
+   read all body data
+   variable amount of info per body, described by ninteger and ndouble
+   to find atoms, must build atom map if not a molecular system
+------------------------------------------------------------------------- */
+
+void ReadData::bodies()
+{
+  int i,m,nchunk,nmax,ninteger,ndouble,tmp,onebody;
+  char *eof;
+
+  int mapflag = 0;
+  if (atom->map_style == 0) {
+    mapflag = 1;
+    atom->map_style = 1;
+    atom->map_init();
+    atom->map_set();
+  }
+
+  // nmax = max # of bodies to read in this chunk
+  // nchunk = actual # readr
+
+  bigint nread = 0;
+  bigint natoms = nbodies;
+
+  while (nread < natoms) {
+    if (natoms-nread > CHUNK) nmax = CHUNK;
+    else nmax = natoms-nread;
+
+    if (me == 0) {
+      nchunk = 0;
+      nlines = 0;
+      m = 0;
+
+      while (nchunk < nmax && nlines <= CHUNK-MAXBODY) {
+        eof = fgets(&buffer[m],MAXLINE,fp);
+        if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
+        sscanf(&buffer[m],"%d %d %d",&tmp,&ninteger,&ndouble);
+        m += strlen(&buffer[m]);
+
+        onebody = 0;
+        if (ninteger) onebody += (ninteger-1)/10 + 1;
+        if (ndouble) onebody += (ndouble-1)/10 + 1;
+        if (onebody+1 > MAXBODY)
+          error->one(FLERR,
+                     "Too many lines in one body in data file - boost MAXBODY");
+
+        for (i = 0; i < onebody; i++) {
+          eof = fgets(&buffer[m],MAXLINE,fp);
+          if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
+          m += strlen(&buffer[m]);
+        }
+
+        nchunk++;
+        nlines += onebody+1;
+      }
+
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+      m++;
+    }
+
+    MPI_Bcast(&nchunk,1,MPI_INT,0,world);
+    MPI_Bcast(&m,1,MPI_INT,0,world);
+    MPI_Bcast(buffer,m,MPI_CHAR,0,world);
+
+    atom->data_bodies(nchunk,buffer,avec_body);
+    nread += nchunk;
+  }
+
+  if (mapflag) {
+    atom->map_delete();
+    atom->map_style = 0;
+  }
+
+  if (me == 0) {
+    if (screen) fprintf(screen,"  " BIGINT_FORMAT " bodies\n",natoms);
+    if (logfile) fprintf(logfile,"  " BIGINT_FORMAT " bodies\n",natoms);
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 void ReadData::bonds()
@@ -1157,6 +1249,7 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
   int ellipsoid_flag = 0;
   int line_flag = 0;
   int tri_flag = 0;
+  int body_flag = 0;
 
   // customize for new sections
   // allocate topology counting vector
@@ -1203,6 +1296,10 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
       if (!avec_tri) error->one(FLERR,"Invalid data file section: Triangles");
       tri_flag = 1;
       skip_lines(ntris);
+    } else if (strcmp(keyword,"Bodies") == 0) {
+      if (!avec_body) error->one(FLERR,"Invalid data file section: Bodies");
+      body_flag = 1;
+      skip_lines(nbodies);
 
     } else if (strcmp(keyword,"Pair Coeffs") == 0) {
       if (force->pair == NULL)
@@ -1434,6 +1531,8 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
     error->one(FLERR,"Needed bonus data not in data file");
   if (ntris && !tri_flag)
     error->one(FLERR,"Needed bonus data not in data file");
+  if (nbodies && !body_flag)
+    error->one(FLERR,"Needed bonus data not in data file");
 }
 
 /* ----------------------------------------------------------------------

src/read_data.h

@@ -49,6 +49,8 @@ class ReadData : protected Pointers {
   class AtomVecLine *avec_line;
   bigint ntris;
   class AtomVecTri *avec_tri;
+  bigint nbodies;
+  class AtomVecBody *avec_body;
 
   void open(char *);
   void scan(int &, int &, int &, int &);
@@ -61,6 +63,7 @@ class ReadData : protected Pointers {
   void atoms();
   void velocities();
   void bonus(bigint, class AtomVec *, const char *);
+  void bodies();
 
   void bonds();
   void angles();

src/read_restart.cpp

@@ -627,6 +627,8 @@ void ReadRestart::header()
       }
 
       atom->create_avec(style,nwords,words);
+      atom->avec->read_restart_settings(fp);
+
       for (int i = 0; i < nwords; i++) delete [] words[i];
       delete [] words;
       delete [] style;

src/write_restart.cpp

@@ -326,6 +326,7 @@ void WriteRestart::header()
   // atom_style must be written before atom class values
   // so read_restart can create class before reading class values
   // if style = hybrid, also write sub-class styles
+  // avec->write_restart() writes atom_style specific info
 
   write_char(ATOM_STYLE,atom->atom_style);
 
@@ -341,6 +342,8 @@ void WriteRestart::header()
     }
   }
 
+  if (me == 0) atom->avec->write_restart_settings(fp);
+
   write_bigint(NATOMS,natoms);
   write_int(NTYPES,atom->ntypes);
   write_bigint(NBONDS,atom->nbonds);