more additions to USER atom styles and bug fixes

src/GRANULAR/fix_pour.cpp

@@ -40,6 +40,7 @@ using namespace MathConst;
 
 enum{ATOM,MOLECULE};
 enum{ONE,RANGE,POLY};
+enum{CONSTANT,EQUAL};    // same as FixGravity
 
 #define EPSILON 0.001
 #define SMALL 1.0e-10
@@ -318,6 +319,10 @@ void FixPour::init()
   if (ifix == -1)
     error->all(FLERR,"No fix gravity defined for fix pour");
 
+  int varflag = ((FixGravity *) modify->fix[ifix])->varflag;
+  if (varflag != CONSTANT) 
+    error->all(FLERR,"Fix gravity for fix pour must be constant");
+
   double xgrav = ((FixGravity *) modify->fix[ifix])->xgrav;
   double ygrav = ((FixGravity *) modify->fix[ifix])->ygrav;
   double zgrav = ((FixGravity *) modify->fix[ifix])->zgrav;

src/MOLECULE/atom_vec_bond.cpp

@@ -60,6 +60,20 @@ AtomVecBond::~AtomVecBond()
   delete [] bond_negative;
 }
 
+/* ----------------------------------------------------------------------
+   grow atom arrays
+   must set local copy of body ptr
+   needed in replicate when 2 atom classes exist and pack_restart() is called
+------------------------------------------------------------------------- */
+
+void AtomVecBond::grow(int n)
+{
+  AtomVec::grow(n);
+  num_bond = atom->num_bond;
+  bond_type = atom->bond_type;
+}
+
+
 /* ----------------------------------------------------------------------
    modify values for AtomVec::pack_restart() to pack
 ------------------------------------------------------------------------- */
@@ -76,8 +90,8 @@ void AtomVecBond::pack_restart_pre(int ilocal)
 
   // flip any negative types to positive and flag which ones
 
-  int *num_bond = atom->num_bond;
+  //int *num_bond = atom->num_bond;
-  int **bond_type = atom->bond_type;
+  //int **bond_type = atom->bond_type;
 
   any_bond_negative = 0;
   for (int m = 0; m < num_bond[ilocal]; m++) {
@@ -98,8 +112,8 @@ void AtomVecBond::pack_restart_post(int ilocal)
   // restore the flagged types to their negative values
 
   if (any_bond_negative) {
-    int *num_bond = atom->num_bond;
+    //int *num_bond = atom->num_bond;
-    int **bond_type = atom->bond_type;
+    //int **bond_type = atom->bond_type;
     for (int m = 0; m < num_bond[ilocal]; m++)
       if (bond_negative[m]) bond_type[ilocal][m] = -bond_type[ilocal][m];
   }

src/MOLECULE/atom_vec_bond.h

@@ -28,6 +28,7 @@ class AtomVecBond : public AtomVec {
  public:
   AtomVecBond(class LAMMPS *);
   ~AtomVecBond();
+  void grow(int);
   void pack_restart_pre(int);
   void pack_restart_post(int);
   void unpack_restart_init(int);
@@ -37,6 +38,9 @@ class AtomVecBond : public AtomVec {
   int any_bond_negative;
   int bond_per_atom;
   int *bond_negative;
+
+  int *num_bond;
+  int **bond_type;
 };
 
 }

src/MOLECULE/atom_vec_full.cpp

@@ -128,7 +128,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
   int *num_improper = atom->num_improper;
   int **improper_type = atom->improper_type;
 
-  int any_bond_negative = 0;
+  any_bond_negative = 0;
   for (int m = 0; m < num_bond[ilocal]; m++) {
     if (bond_type[ilocal][m] < 0) {
       bond_negative[m] = 1;
@@ -137,7 +137,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
     } else bond_negative[m] = 0;
   }
 
-  int any_angle_negative = 0;
+  any_angle_negative = 0;
   for (int m = 0; m < num_angle[ilocal]; m++) {
     if (angle_type[ilocal][m] < 0) {
       angle_negative[m] = 1;
@@ -146,7 +146,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
     } else angle_negative[m] = 0;
   }
 
-  int any_dihedral_negative = 0;
+  any_dihedral_negative = 0;
   for (int m = 0; m < num_dihedral[ilocal]; m++) {
     if (dihedral_type[ilocal][m] < 0) {
       dihedral_negative[m] = 1;
@@ -155,7 +155,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
     } else dihedral_negative[m] = 0;
   }
 
-  int any_improper_negative = 0;
+  any_improper_negative = 0;
   for (int m = 0; m < num_improper[ilocal]; m++) {
     if (improper_type[ilocal][m] < 0) {
       improper_negative[m] = 1;

src/MOLECULE/atom_vec_molecular.cpp

@@ -128,7 +128,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
   int *num_improper = atom->num_improper;
   int **improper_type = atom->improper_type;
 
-  int any_bond_negative = 0;
+  any_bond_negative = 0;
   for (int m = 0; m < num_bond[ilocal]; m++) {
     if (bond_type[ilocal][m] < 0) {
       bond_negative[m] = 1;
@@ -137,7 +137,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
     } else bond_negative[m] = 0;
   }
 
-  int any_angle_negative = 0;
+  any_angle_negative = 0;
   for (int m = 0; m < num_angle[ilocal]; m++) {
     if (angle_type[ilocal][m] < 0) {
       angle_negative[m] = 1;
@@ -146,7 +146,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
     } else angle_negative[m] = 0;
   }
 
-  int any_dihedral_negative = 0;
+  any_dihedral_negative = 0;
   for (int m = 0; m < num_dihedral[ilocal]; m++) {
     if (dihedral_type[ilocal][m] < 0) {
       dihedral_negative[m] = 1;
@@ -155,7 +155,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
     } else dihedral_negative[m] = 0;
   }
 
-  int any_improper_negative = 0;
+  any_improper_negative = 0;
   for (int m = 0; m < num_improper[ilocal]; m++) {
     if (improper_type[ilocal][m] < 0) {
       improper_negative[m] = 1;

src/USER-AWPMD/atom_vec_wavepacket.h

@@ -11,11 +11,6 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-/* ----------------------------------------------------------------------
-   Contributing author: Ilya Valuev (JIHT RAS)
-------------------------------------------------------------------------- */
-
-
 #ifdef ATOM_CLASS
 
 AtomStyle(wavepacket,AtomVecWavepacket)
@@ -32,77 +27,11 @@ namespace LAMMPS_NS {
 class AtomVecWavepacket : public AtomVec {
 public:
   AtomVecWavepacket(class LAMMPS *);
-  ~AtomVecWavepacket() {}
-  void grow(int);
-  void grow_reset();
-  void copy(int, int, int);
   void force_clear(int, size_t);
-  int pack_comm(int, int *, double *, int, int *);
+  void create_atom_post(int);
-  int pack_comm_vel(int, int *, double *, int, int *);
+  void data_atom_post(int);
-  int pack_comm_hybrid(int, int *, double *);
-  void unpack_comm(int, int, double *);
-  void unpack_comm_vel(int, int, double *);
-  int unpack_comm_hybrid(int, int, double *);
-  int pack_reverse(int, int, double *);
-  int pack_reverse_hybrid(int, int, double *);
-  void unpack_reverse(int, int *, double *);
-  int unpack_reverse_hybrid(int, int *, double *);
-  int pack_border(int, int *, double *, int, int *);
-  int pack_border_vel(int, int *, double *, int, int *);
-  int pack_border_hybrid(int, int *, double *);
-  void unpack_border(int, int, double *);
-  void unpack_border_vel(int, int, double *);
-  int unpack_border_hybrid(int, int, double *);
-  int pack_exchange(int, double *);
-  int unpack_exchange(double *);
-  int size_restart();
-  int pack_restart(int, double *);
-  int unpack_restart(double *);
-  void create_atom(int, double *);
-  void data_atom(double *, imageint, char **);
-  int data_atom_hybrid(int, char **);
-  void data_vel(int, char **);
-  int data_vel_hybrid(int, char **);
-  void pack_data(double **);
-  int pack_data_hybrid(int, double *);
-  void write_data(FILE *, int, double **);
-  int write_data_hybrid(FILE *, double *);
-  void pack_vel(double **);
-  int pack_vel_hybrid(int, double *);
-  void write_vel(FILE *, int, double **);
-  int write_vel_hybrid(FILE *, double *);
   int property_atom(char *);
   void pack_property_atom(int, double *, int, int);
-  bigint memory_usage();
-
-private:
-  tagint *tag;
-  int *type,*mask;
-  imageint *image;
-  double **x,**v,**f;
-
-  ///\en spin: -1 or 1 for electron, 0 for ion (compatible with eff)
-  int *spin;
-  ///\en charge: must be specified in the corresponding units (-1 for electron in real units, eff compatible)
-  double *q;
-  ///\en width of the wavepacket (compatible with eff)
-  double *eradius;
-  ///\en width velocity for the wavepacket (compatible with eff)
-  double *ervel;
-  ///\en (generalized) force on width  (compatible with eff)
-  double *erforce;
-
-  // AWPMD- specific:
-  ///\en electron tag: must be the same for the WPs belonging to the same electron
-  int *etag;
-  ///\en wavepacket split coefficients: cre, cim, size is 2*N
-  double *cs;
-  ///\en force on wavepacket split coefficients: re, im, size is 2*N
-  double *csforce;
-  ///\en (generalized) force on velocity, size is 3*N
-  double *vforce;
-   ///\en (generalized) force on radius velocity, size is N
-  double *ervelforce;
 };
 
 }

src/USER-AWPMD/fix_nve_awpmd.cpp

@@ -74,8 +74,6 @@ void FixNVEAwpmd::init()
 
 void FixNVEAwpmd::initial_integrate(int /* vflag */)
 {
-
-
   // update v,vr and x,radius of atoms in group
 
   double **x = atom->x;
@@ -84,7 +82,7 @@ void FixNVEAwpmd::initial_integrate(int /* vflag */)
   double *ervel = atom->ervel;
   double **f = atom->f;
   double *erforce = atom->erforce;
-  double *vforce=atom->vforce;
+  double **vforce=atom->vforce;
   double *ervelforce=atom->ervelforce;
 
   double *mass = atom->mass;
@@ -101,7 +99,7 @@ void FixNVEAwpmd::initial_integrate(int /* vflag */)
       double dtfm = dtf / mass[type[i]];
       double dtfmr=dtfm;
       for(int j=0;j<3;j++){
-        x[i][j] += dtv*vforce[3*i+j];
+        x[i][j] += dtv*vforce[i][j];
         v[i][j] += dtfm*f[i][j];
       }
       eradius[i]+= dtv*ervelforce[i];

src/USER-DPD/atom_vec_dpd.cpp

@@ -47,7 +47,7 @@ AtomVecDPD::AtomVecDPD(LAMMPS *lmp) : AtomVec(lmp)
   fields_restart = (char *) "dpdTheta uCond uMech uChem";
   fields_create = (char *) "rho dpdTheta uCond uMech uChem uCG uCGnew duChem";
   fields_data_atom = (char *) "id type dpdTheta x";
-  fields_data_vel = (char *) "id v omega";
+  fields_data_vel = (char *) "id v";
 
   setup_fields();
 }

src/USER-EFF/atom_vec_electron.cpp

@@ -44,8 +44,6 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp)
   molecular = 0;
   forceclearflag = 1;
 
-  atom->ecp_flag = 0;
-
   atom->electron_flag = 1;
   atom->q_flag = atom->spin_flag = atom->eradius_flag =
     atom->ervel_flag = atom->erforce_flag = 1;
@@ -99,7 +97,6 @@ void AtomVecElectron::create_atom_post(int ilocal)
 void AtomVecElectron::data_atom_post(int ilocal)
 {
   atom->ervel[ilocal] = 0.0;
-  if (atom->spin[ilocal] == 3) atom->ecp_flag = 1;
 }
 
 /* ----------------------------------------------------------------------

src/USER-EFF/pair_eff_cut.cpp

@@ -801,7 +801,7 @@ void PairEffCut::settings(int narg, char **arg)
 
   int atype;
   int iarg = 1;
-  int ecp_found = 0;
+  ecp_found = 0;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"limit/eradius") == 0) {
@@ -821,17 +821,15 @@ void PairEffCut::settings(int narg, char **arg)
         else if (strcmp(arg[iarg+1],"O") == 0) ecp_type[atype] = 8;
         else if (strcmp(arg[iarg+1],"Al") == 0) ecp_type[atype] = 13;
         else if (strcmp(arg[iarg+1],"Si") == 0) ecp_type[atype] = 14;
-        else error->all(FLERR, "Note: there are no default parameters for this atom ECP\n");
+        else error->all(FLERR, "No default parameters for this atom ECP\n");
         iarg += 2;
         ecp_found = 1;
       }
-    }
+    } else error->all(FLERR,"Illegal pair style command");
   }
 
-  if (!ecp_found && atom->ecp_flag)
-    error->all(FLERR,"Need to specify ECP type on pair_style command");
-
   // Need to introduce 2 new constants w/out changing update.cpp
+
   if (force->qqr2e==332.06371) {        // i.e. Real units chosen
     h2e = 627.509;                      // hartree->kcal/mol
     hhmss2e = 175.72044219620075;       // hartree->kcal/mol * (Bohr->Angstrom)^2
@@ -872,9 +870,24 @@ void PairEffCut::init_style()
 
   if (update->whichflag == 1) {
     if (force->qqr2e == 332.06371 && update->dt == 1.0)
-      error->all(FLERR,"You must lower the default real units timestep for pEFF ");
+      error->all(FLERR,"Must lower the default real units timestep for pEFF ");
   }
 
+  // check if any atom's spin = 3 and ECP type was not set
+
+  int *spin = atom->spin;
+  int nlocal = atom->nlocal;
+
+  int flag = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (spin[i] == 3) flag = 1;
+
+  int flagall;
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+
+  if (flagall && !ecp_found)
+    error->all(FLERR,"Need to specify ECP type on pair_style command");
+
   // need a half neigh list and optionally a granular history neigh list
 
   neighbor->request(this,instance_me);

src/USER-EFF/pair_eff_cut.h

@@ -46,10 +46,12 @@ class PairEffCut : public Pair {
 
  private:
   int limit_eradius_flag, pressure_with_evirials_flag;
+  int ecp_found;
   double cut_global;
   double **cut;
   int ecp_type[100];
-  double PAULI_CORE_A[100], PAULI_CORE_B[100], PAULI_CORE_C[100], PAULI_CORE_D[100], PAULI_CORE_E[100];
+  double PAULI_CORE_A[100],PAULI_CORE_B[100],PAULI_CORE_C[100];
+  double PAULI_CORE_D[100],PAULI_CORE_E[100];
   double hhmss2e, h2e;
 
   int nmax;

src/USER-MESO/atom_vec_edpd.cpp

@@ -41,16 +41,16 @@ AtomVecEDPD::AtomVecEDPD(LAMMPS *lmp) : AtomVec(lmp)
   // order of fields in a string does not matter
   // except: fields_data_atom & fields_data_vel must match data file
 
-  fields_grow = (char *) "edpd_cv edpd_temp edpd_flux vest";
+  fields_grow = (char *) "edpd_cv edpd_temp edpd_flux vest vest_temp";
-  fields_copy = (char *) "edpd_cv edpd_temp edpd_flux vest";
+  fields_copy = (char *) "edpd_cv edpd_temp edpd_flux vest vest_temp";
-  fields_comm = (char *) "edpd_temp vest";
+  fields_comm = (char *) "edpd_temp vest vest_temp";
-  fields_comm_vel = (char *) "edpd_temp vest";
+  fields_comm_vel = (char *) "edpd_temp vest vest_temp";
   fields_reverse = (char *) "edpd_flux";
-  fields_border = (char *) "edpd_cv edpd_temp vest";
+  fields_border = (char *) "edpd_cv edpd_temp vest vest_temp";
-  fields_border_vel = (char *) "edpd_cv edpd_temp vest";
+  fields_border_vel = (char *) "edpd_cv edpd_temp vest vest_temp";
-  fields_exchange = (char *) "edpd_cv edpd_temp vest";
+  fields_exchange = (char *) "edpd_cv edpd_temp vest vest_temp";
-  fields_restart = (char * ) "edpd_cv edpd_temp vest";
+  fields_restart = (char * ) "edpd_cv edpd_temp vest vest_temp";
-  fields_create = (char *) "edpd_cv edpd_temp edpd_flux vest";
+  fields_create = (char *) "edpd_cv edpd_temp edpd_flux vest vest_temp";
   fields_data_atom = (char *) "id type edpd_temp edpd_cv x";
   fields_data_vel = (char *) "id v";
 
@@ -85,7 +85,7 @@ void AtomVecEDPD::create_atom_post(int ilocal)
 {
   atom->edpd_temp[ilocal] = 1.0;
   atom->edpd_cv[ilocal]= 1.0e5;
-  atom->vest[ilocal][3] = atom->edpd_temp[ilocal];
+  atom->vest_temp[ilocal] = atom->edpd_temp[ilocal];
 }
 
 /* ----------------------------------------------------------------------
@@ -99,5 +99,5 @@ void AtomVecEDPD::data_atom_post(int ilocal)
   atom->vest[ilocal][0] = 0.0;
   atom->vest[ilocal][1] = 0.0;
   atom->vest[ilocal][2] = 0.0;
-  atom->vest[ilocal][3] = atom->edpd_temp[ilocal];
+  atom->vest_temp[ilocal] = atom->edpd_temp[ilocal];
 }

src/USER-MESO/atom_vec_tdpd.cpp

@@ -63,7 +63,7 @@ void AtomVecTDPD::process_args(int narg, char **arg)
   cc_species = atom->cc_species;
 
   atom->add_peratom_change_columns("cc",cc_species);
-  atom->add_peratom_change_columns("cc_species",cc_species);
+  atom->add_peratom_change_columns("cc_flux",cc_species);
 
   // delay setting up of fields until now
 

src/USER-MESO/fix_mvv_edpd.cpp

@@ -88,6 +88,7 @@ void FixMvvEDPD::initial_integrate(int /*vflag*/)
   double *edpd_flux = atom->edpd_flux;
   double *edpd_cv = atom->edpd_cv;
   double **vest = atom->vest;
+  double *vest_temp = atom->vest_temp;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
@@ -105,7 +106,7 @@ void FixMvvEDPD::initial_integrate(int /*vflag*/)
      vest[i][0] = v[i][0] + dtfm * f[i][0];
      vest[i][1] = v[i][1] + dtfm * f[i][1];
      vest[i][2] = v[i][2] + dtfm * f[i][2];
-     vest[i][3] = edpd_temp[i] + dtT * edpd_flux[i];
+     vest_temp[i] = edpd_temp[i] + dtT * edpd_flux[i];
 
      x[i][0] += dtv * vest[i][0];
      x[i][1] += dtv * vest[i][1];
@@ -131,6 +132,7 @@ void FixMvvEDPD::final_integrate()
   double *edpd_flux = atom->edpd_flux;
   double *edpd_cv = atom->edpd_cv;
   double **vest = atom->vest;
+  double *vest_temp = atom->vest_temp;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
@@ -148,7 +150,7 @@ void FixMvvEDPD::final_integrate()
      v[i][0] = vest[i][0] + dtfm * f[i][0];
      v[i][1] = vest[i][1] + dtfm * f[i][1];
      v[i][2] = vest[i][2] + dtfm * f[i][2];
-     edpd_temp[i] = vest[i][3] + dtT * edpd_flux[i];
+     edpd_temp[i] = vest_temp[i] + dtT * edpd_flux[i];
   }
 }
 

src/USER-SMD/atom_vec_smd.cpp

@@ -63,7 +63,7 @@ AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp)
 
   fields_grow = (char *) 
     "de vfrac rmass x0 radius contact_radius molecule "
-    "smd_data_9 e vest tlsph_stress "
+    "smd_data_9 e vest smd_stress "
     "eff_plastic_strain eff_plastic_strain_rate damage";
   fields_copy = (char *) 
     "vfrac rmass x0 radius contact_radius molecule e "
@@ -91,7 +91,7 @@ AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp)
     "eff_plastic_strain eff_plastic_strain_rate smd_data_9 smd_stress damage";
   fields_data_atom = (char *) 
     "id type molecule vfrac rmass radius contact_radius x";
-  fields_data_vel = (char *) "id v vest";
+  fields_data_vel = (char *) "id v";
 
   // set these array sizes based on defines
 

src/atom.cpp

@@ -566,6 +566,7 @@ void Atom::peratom_create()
 
   add_peratom("edpd_cv",&edpd_cv,DOUBLE,0);
   add_peratom("edpd_temp",&edpd_temp,DOUBLE,0);
+  add_peratom("vest_temp",&vest_temp,DOUBLE,0);
   add_peratom("edpd_flux",&edpd_flux,DOUBLE,0,1);     // set per-thread flag
   add_peratom("cc",&cc,DOUBLE,1);
   add_peratom("cc_flux",&cc_flux,DOUBLE,1,1);         // set per-thread flag
@@ -625,8 +626,11 @@ void Atom::add_peratom(const char *name, void *address,
 
 void Atom::add_peratom_change_columns(const char *name, int cols)
 {
+  int i;
   for (int i = 0; i < nperatom; i++)
     if (strcmp(name,peratom[i].name) == 0) peratom[i].cols = cols;
+  if (i == nperatom) 
+    error->all(FLERR,"Could not find name of peratom array for column change");
 }
 
 /* ----------------------------------------------------------------------

src/atom.h

@@ -124,6 +124,7 @@ class Atom : protected Pointers {
 
   double **cc,**cc_flux;         // cc = chemical concentration
   double *edpd_temp,*edpd_flux;  // temperature and heat flux
+  double *vest_temp;
   double *edpd_cv;               // heat capacity
   int cc_species;
 
@@ -162,10 +163,6 @@ class Atom : protected Pointers {
 
   int sp_flag;
 
-  // USER-EFF package
-
-  int ecp_flag;
-
   // USER-SMD package
 
   int smd_flag;

src/atom_vec.cpp

@@ -45,6 +45,7 @@ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp)
   forceclearflag = 0;
   maxexchange = 0;
   bonus_flag = 0;
+  size_forward_bonus = size_border_bonus = 0;
 
   kokkosable = 0;
 
@@ -176,7 +177,7 @@ int AtomVec::grow_nmax_bonus(int nmax_bonus)
 
 void AtomVec::grow(int n)
 {
-  int i,datatype,cols,maxcols;
+  int datatype,cols,maxcols;
   void *pdata,*plength;
 
   if (n == 0) grow_nmax();
@@ -193,7 +194,7 @@ void AtomVec::grow(int n)
   v = memory->grow(atom->v,nmax,3,"atom:v");
   f = memory->grow(atom->f,nmax*nthreads,3,"atom:f");
 
-  for (i = 0; i < ngrow; i++) {
+  for (int i = 0; i < ngrow; i++) {
     pdata = mgrow.pdata[i];
     datatype = mgrow.datatype[i];
     cols = mgrow.cols[i];
@@ -415,7 +416,7 @@ int AtomVec::pack_comm(int n, int *list, double *buf,
 /* ---------------------------------------------------------------------- */
 
 int AtomVec::pack_comm_vel(int n, int *list, double *buf,
-                                 int pbc_flag, int *pbc)
+                           int pbc_flag, int *pbc)
 {
   int i,j,m,mm,nn,datatype,cols;
   double dx,dy,dz,dvx,dvy,dvz;
@@ -1375,15 +1376,15 @@ int AtomVec::size_restart()
 
   if (nrestart) {
     for (nn = 0; nn < nrestart; nn++) {
-      cols = mrestart.cols[i];
+      cols = mrestart.cols[nn];
       if (cols == 0) n += nlocal;
       else if (cols > 0) n += cols*nlocal;
       else {
         collength = mrestart.collength[nn];
         plength = mrestart.plength[nn];
         for (i = 0; i < nlocal; i++) {
-          if (collength) ncols = (*((int ***) plength))[nlocal][collength-1];
+          if (collength) ncols = (*((int ***) plength))[i][collength-1];
-          else ncols = (*((int **) plength))[nlocal];
+          else ncols = (*((int **) plength))[i];
           n += ncols;
         }
       }
@@ -1435,7 +1436,7 @@ int AtomVec::pack_restart(int i, double *buf)
       if (cols == 0) {
         double *vec = *((double **) pdata);
         buf[m++] = vec[i];
-      } else if (ncols > 0) {
+      } else if (cols > 0) {
         double **array = *((double ***) pdata);
         for (mm = 0; mm < cols; mm++)
           buf[m++] = array[i][mm];
@@ -1469,7 +1470,7 @@ int AtomVec::pack_restart(int i, double *buf)
       if (cols == 0) {
         bigint *vec = *((bigint **) pdata);
         buf[m++] = ubuf(vec[i]).d;
-      } else if (ncols > 0) {
+      } else if (cols > 0) {
         bigint **array = *((bigint ***) pdata);
         for (mm = 0; mm < cols; mm++)
           buf[m++] = ubuf(array[i][mm]).d;
@@ -1705,7 +1706,7 @@ void AtomVec::data_atom(double *coord, imageint imagetmp, char **values)
         vec[nlocal] = utils::numeric(FLERR,values[ivalue++],true,lmp);
       } else {
         double **array = *((double ***) pdata);
-        if (array == atom->x) {      // already set by coord arg
+        if (array == atom->x) {      // x was already set by coord arg
           ivalue += cols;
           continue;
         }
@@ -1878,9 +1879,9 @@ void AtomVec::data_vel(int ilocal, char **values)
   v[ilocal][1] = utils::numeric(FLERR,values[1],true,lmp);
   v[ilocal][2] = utils::numeric(FLERR,values[2],true,lmp);
 
-  if (ndata_vel) {
+  if (ndata_vel > 2) {
     int ivalue = 3;
-    for (n = 0; n < ndata_vel; n++) {
+    for (n = 2; n < ndata_vel; n++) {
       pdata = mdata_vel.pdata[n];
       datatype = mdata_vel.datatype[n];
       cols = mdata_vel.cols[n];
@@ -2434,6 +2435,7 @@ void AtomVec::setup_fields()
     else size_data_atom += cols;
   }
 
+
   size_data_vel = 0;
   for (n = 0; n < ndata_vel; n++) {
     cols = mdata_vel.cols[n];

src/atom_vec.h

@@ -74,7 +74,7 @@ class AtomVec : protected Pointers {
 
   virtual void force_clear(int, size_t) {}
 
-  void grow(int);
+  virtual void grow(int);
   void copy(int, int, int);
 
   virtual void copy_bonus(int, int, int) {}
@@ -122,13 +122,11 @@ class AtomVec : protected Pointers {
 
   void data_atom(double *, imageint, char **);
   virtual void data_atom_post(int) {}
-
+  virtual void data_atom_bonus(int, char **) {}
-  void data_atom_bonus(int, char **) {}
+  virtual void data_body(int, int, int, int *, double *) {}
-  void data_body(int, int, int, int *, double *) {}
 
   void pack_data(double **);
   void write_data(FILE *, int, double **);
-
   virtual void pack_data_pre(int) {}
   virtual void pack_data_post(int) {}
 
@@ -145,19 +143,19 @@ class AtomVec : protected Pointers {
   int pack_improper(tagint **);
   void write_improper(FILE *, int, tagint **, int);
 
-  int property_atom(char *) {return -1;}
+  virtual int property_atom(char *) {return -1;}
-  void pack_property_atom(int, double *, int, int) {}
+  virtual void pack_property_atom(int, double *, int, int) {}
 
   bigint memory_usage();
   virtual bigint memory_usage_bonus() {}
 
  protected:
-  int nmax;                             // local copy of atom->nmax
+  int nmax;                    // local copy of atom->nmax
-  int deform_vremap;                    // local copy of domain properties
+  int deform_vremap;           // local copy of domain properties
   int deform_groupbit;
   double *h_rate;
 
-  tagint *tag;                          // peratom fields common to all styles
+  tagint *tag;                 // peratom fields common to all styles
   int *type,*mask;
   imageint *image;
   double **x,**v,**f;

src/atom_vec_body.cpp

@@ -40,8 +40,8 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
   // size_data_bonus is not used by Atom class for body style
 
   size_forward_bonus = 4;
-  size_border_bonus = 9;
+  size_border_bonus = 10;
-  size_restart_bonus_one = 9;
+  size_restart_bonus_one = 10;
   size_data_bonus = 0;
 
   atom->body_flag = 1;
@@ -72,11 +72,9 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
   fields_border_vel = (char *) "radius rmass angmom";
   fields_exchange = (char *) "radius rmass angmom";
   fields_restart = (char *) "radius rmass angmom";
-  fields_create = (char *) "radius rmass angmom tri";
+  fields_create = (char *) "radius rmass angmom body";
   fields_data_atom = (char *) "id type body rmass x";
   fields_data_vel = (char *) "id v angmom";
-
-  setup_fields();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -127,6 +125,20 @@ void AtomVecBody::process_args(int narg, char **arg)
 
   size_forward_bonus += bptr->size_forward;
   size_border_bonus += bptr->size_border;
+
+  setup_fields();
+}
+
+/* ----------------------------------------------------------------------
+   grow atom arrays
+   must set local copy of body ptr
+   needed in replicate when 2 atom classes exist and pack_restart() is called
+------------------------------------------------------------------------- */
+
+void AtomVecBody::grow(int n)
+{
+  AtomVec::grow(n);
+  body = atom->body;
 }
 
 /* ----------------------------------------------------------------------
@@ -150,8 +162,6 @@ void AtomVecBody::grow_bonus()
 
 void AtomVecBody::copy_bonus(int i, int j, int delflag)
 {
-  int *body = atom->body;
-
   // if deleting atom J via delflag and J has bonus data, then delete it
 
   if (delflag && body[j] >= 0) {
@@ -206,8 +216,6 @@ int AtomVecBody::pack_comm_bonus(int n, int *list, double *buf)
   int i,j,m;
   double *quat;
 
-  int *body = atom->body;
-
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
@@ -231,8 +239,6 @@ void AtomVecBody::unpack_comm_bonus(int n, int first, double *buf)
   int i,m,last;
   double *quat;
 
-  int *body = atom->body;
-
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
@@ -254,8 +260,6 @@ int AtomVecBody::pack_border_bonus(int n, int *list, double *buf)
   int i,j,m;
   double *quat,*inertia;
 
-  int *body = atom->body;
-
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
@@ -287,8 +291,6 @@ int AtomVecBody::unpack_border_bonus(int n, int first, double *buf)
   int i,j,m,last;
   double *quat,*inertia;
 
-  int *body = atom->body;
-
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
@@ -330,8 +332,6 @@ int AtomVecBody::pack_exchange_bonus(int i, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   if (body[i] < 0) buf[m++] = ubuf(0).d;
   else {
     buf[m++] = ubuf(1).d;
@@ -363,8 +363,6 @@ int AtomVecBody::unpack_exchange_bonus(int ilocal, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   body[ilocal] = (int) ubuf(buf[m++]).i;
   if (body[ilocal] == 0) body[ilocal] = -1;
   else {
@@ -408,8 +406,6 @@ int AtomVecBody::size_restart_bonus()
 {
   int i;
 
-  int *body = atom->body;
-
   int n = 0;
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
@@ -418,8 +414,7 @@ int AtomVecBody::size_restart_bonus()
       if (intdoubleratio == 1) n += bonus[body[i]].ninteger;
       else n += (bonus[body[i]].ninteger+1)/2;
       n += bonus[body[i]].ndouble;
-    }
+    } else n++;
-    n++;
   }
 
   return n;
@@ -435,8 +430,6 @@ int AtomVecBody::pack_restart_bonus(int i, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   if (body[i] < 0) buf[m++] = ubuf(0).d;
   else {
     buf[m++] = ubuf(1).d;
@@ -470,8 +463,6 @@ int AtomVecBody::unpack_restart_bonus(int ilocal, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   body[ilocal] = (int) ubuf(buf[m++]).i;
   if (body[ilocal] == 0) body[ilocal] = -1;
   else {

src/atom_vec_body.h

@@ -43,6 +43,7 @@ class AtomVecBody : public AtomVec {
   ~AtomVecBody();
   void process_args(int, char **);
 
+  void grow(int);
   void copy_bonus(int, int, int);
   void clear_bonus();
   int pack_comm_bonus(int, int *, double *);
@@ -74,6 +75,8 @@ class AtomVecBody : public AtomVec {
   int intdoubleratio;            // sizeof(double) / sizeof(int)
   int body_flag;
 
+  int *body;
+
   MyPoolChunk<int> *icp;
   MyPoolChunk<double> *dcp;
 

src/atom_vec_ellipsoid.cpp

@@ -38,7 +38,7 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp)
 
   size_forward_bonus = 4;
   size_border_bonus = 8;
-  size_restart_bonus_one = 7;
+  size_restart_bonus_one = 8;
   size_data_bonus = 8;
 
   atom->ellipsoid_flag = 1;
@@ -320,7 +320,7 @@ int AtomVecEllipsoid::size_restart_bonus()
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
     if (ellipsoid[i] >= 0) n += size_restart_bonus_one;
-    n++;
+    else n++;
   }
 
   return n;

src/atom_vec_line.cpp

@@ -37,7 +37,7 @@ AtomVecLine::AtomVecLine(LAMMPS *lmp) : AtomVec(lmp)
 
   size_forward_bonus = 1;
   size_border_bonus = 3;
-  size_restart_bonus_one = 2;
+  size_restart_bonus_one = 3;
   size_data_bonus = 5;
 
   atom->line_flag = 1;
@@ -286,7 +286,7 @@ int AtomVecLine::size_restart_bonus()
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
     if (line[i] >= 0) n += size_restart_bonus_one;
-    n++;
+    else n++;
   }
 
   return n;

src/atom_vec_sphere.cpp

@@ -61,18 +61,22 @@ AtomVecSphere::AtomVecSphere(LAMMPS *lmp) : AtomVec(lmp)
 
 void AtomVecSphere::process_args(int narg, char **arg)
 {
-  if (narg == 0) return;
+  if (narg != 0 && narg != 1) 
-  if (narg != 1) error->all(FLERR,"Illegal atom_style sphere command");
-
-  radvary = utils::numeric(FLERR,arg[0],true,lmp);
-  if (radvary < 0 || radvary > 1)
     error->all(FLERR,"Illegal atom_style sphere command");
-  if (radvary == 0) return;
+
+  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 command");
+  }
 
   // dynamic particle radius and mass must be communicated every step
 
-  fields_comm = (char *) "radius rmass";
+  if (radvary) {
-  fields_comm_vel = (char *) "radius rmass omega";
+    fields_comm = (char *) "radius rmass";
+    fields_comm_vel = (char *) "radius rmass omega";
+  }
 
   // delay setting up of fields until now
 

src/atom_vec_tri.cpp

@@ -38,7 +38,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
 
   size_forward_bonus = 4;
   size_border_bonus = 17;
-  size_restart_bonus_one = 16;
+  size_restart_bonus_one = 17;
   size_data_bonus = 10;
 
   atom->tri_flag = 1;
@@ -64,7 +64,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
   fields_border_vel = (char *) "molecule radius rmass omega";
   fields_exchange = (char *) "molecule radius rmass omega angmom";
   fields_restart = (char *) "molecule radius rmass omega angmom";
-  fields_create = (char *) "molecule radius rmass omega angmom line";
+  fields_create = (char *) "molecule radius rmass omega angmom tri";
   fields_data_atom = (char *) "id molecule type tri rmass x";
   fields_data_vel = (char *) "id v omega angmom";
 
@@ -381,7 +381,7 @@ int AtomVecTri::size_restart_bonus()
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
     if (tri[i] >= 0) n += size_restart_bonus_one;
-    n++;
+    else n++;
   }
 
   return n;

src/fix_gravity.cpp

@@ -31,7 +31,7 @@ using namespace FixConst;
 using namespace MathConst;
 
 enum{CHUTE,SPHERICAL,VECTOR};
-enum{CONSTANT,EQUAL};
+enum{CONSTANT,EQUAL};          // same as FixPour
 
 /* ---------------------------------------------------------------------- */
 
@@ -134,6 +134,16 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
 
   eflag = 0;
   egrav = 0.0;
+
+  // set gravity components once and for all if CONSTANT
+
+  varflag = CONSTANT;
+  if (mstyle != CONSTANT || vstyle != CONSTANT || pstyle != CONSTANT ||
+      tstyle != CONSTANT || xstyle != CONSTANT || ystyle != CONSTANT ||
+      zstyle != CONSTANT) varflag = EQUAL;
+
+  if (varflag == CONSTANT) set_acceleration();
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -222,15 +232,6 @@ void FixGravity::init()
     if (!input->variable->equalstyle(zvar))
       error->all(FLERR,"Variable for fix gravity is invalid style");
   }
-
-  varflag = CONSTANT;
-  if (mstyle != CONSTANT || vstyle != CONSTANT || pstyle != CONSTANT ||
-      tstyle != CONSTANT || xstyle != CONSTANT || ystyle != CONSTANT ||
-      zstyle != CONSTANT) varflag = EQUAL;
-
-  // set gravity components once and for all
-
-  if (varflag == CONSTANT) set_acceleration();
 }
 
 /* ---------------------------------------------------------------------- */

src/replicate.cpp

@@ -331,35 +331,43 @@ void Replicate::command(int narg, char **arg)
 
     if (me == 0 && screen) {
       fprintf(screen,"  bounding box image = (%i %i %i) to (%i %i %i)\n",
-              _imagelo[0],_imagelo[1],_imagelo[2],_imagehi[0],_imagehi[1],_imagehi[2]);
+              _imagelo[0],_imagelo[1],_imagelo[2],
+              _imagehi[0],_imagehi[1],_imagehi[2]);
       fprintf(screen,"  bounding box extra memory = %.2f MB\n",
               (double)size_buf_all*sizeof(double)/1024/1024);
     }
 
     // rnk offsets
 
-    int * disp_buf_rnk;
+    int *disp_buf_rnk;
     memory->create(disp_buf_rnk, nprocs, "replicate:disp_buf_rnk");
     disp_buf_rnk[0] = 0;
-    for (i=1; i<nprocs; ++i) disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1];
+    for (i = 1; i < nprocs; i++) 
+      disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1];
 
     // allgather buf_all
 
     double * buf_all;
     memory->create(buf_all, size_buf_all, "replicate:buf_all");
 
-    MPI_Allgatherv(buf, n, MPI_DOUBLE, buf_all, size_buf_rnk, disp_buf_rnk, MPI_DOUBLE, world);
+    MPI_Allgatherv(buf,n,MPI_DOUBLE,buf_all,size_buf_rnk,disp_buf_rnk,
+                   MPI_DOUBLE,world);
 
     // bounding box of original unwrapped system
 
     double _orig_lo[3], _orig_hi[3];
     if (triclinic) {
-      _orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd + _imagelo[1] * old_xy + _imagelo[2] * old_xz;
+      _orig_lo[0] = domain->boxlo[0] + 
-      _orig_lo[1] = domain->boxlo[1] + _imagelo[1] * old_yprd + _imagelo[2] * old_yz;
+        _imagelo[0] * old_xprd + _imagelo[1] * old_xy + _imagelo[2] * old_xz;
+      _orig_lo[1] = domain->boxlo[1] + 
+        _imagelo[1] * old_yprd + _imagelo[2] * old_yz;
       _orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd;
 
-      _orig_hi[0] = domain->boxlo[0] + (_imagehi[0]+1) * old_xprd + (_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz;
+      _orig_hi[0] = domain->boxlo[0] + 
-      _orig_hi[1] = domain->boxlo[1] + (_imagehi[1]+1) * old_yprd + (_imagehi[2]+1) * old_yz;
+        (_imagehi[0]+1) * old_xprd + 
+        (_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz;
+      _orig_hi[1] = domain->boxlo[1] + 
+        (_imagehi[1]+1) * old_yprd + (_imagehi[2]+1) * old_yz;
       _orig_hi[2] = domain->boxlo[2] + (_imagehi[2]+1) * old_zprd;
     } else {
       _orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd;
@@ -605,7 +613,8 @@ void Replicate::command(int narg, char **arg)
     MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world);
     double avg = (double) sum / nprocs;
     if (me == 0 && screen)
-      fprintf(screen,"  average # of replicas added to proc = %.2f out of %i (%.2f %%)\n",
+      fprintf(screen,"  average # of replicas added to proc = %.2f "
+              "out of %i (%.2f %%)\n",
               avg,nx*ny*nz,avg/(nx*ny*nz)*100.0);
 
   } else {

src/verlet.cpp

@@ -109,7 +109,7 @@ void Verlet::setup(int flag)
   domain->pbc();
   domain->reset_box();
   comm->setup();
-  if (neighbor->style) neighbor->setup_bins();
+  if (neighbor->style) neighbor->setup_bins(); 
   comm->exchange();
   if (atom->sortfreq > 0) atom->sort();
   comm->borders();