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

src/ASPHERE/atom_vec_ellipsoid.cpp

@@ -37,19 +37,21 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) :
   AtomVec(lmp, narg, arg)
 {
   molecular = 0;
-  mass_type = 1;
-  shape_type = 1;
 
   comm_x_only = comm_f_only = 0;
   size_forward = 7;
   size_reverse = 6;
-  size_border = 10;
+  size_border = 13;
   size_velocity = 6;
-  size_data_atom = 9;
+  size_data_atom = 13;
   size_data_vel = 7;
-  xcol_data = 3;
+  xcol_data = 7;
 
-  atom->angmom_flag = atom->torque_flag = atom->quat_flag = 1;
+  atom->ellipsoid_flag = 1;
+  atom->shape_flag = atom->rmass_flag = atom->quat_flag = 
+    atom->angmom_flag = atom->torque_flag = 1;
+
+  PI = 4.0*atan(1.0);
 }
 
 /* ----------------------------------------------------------------------
@@ -74,6 +76,8 @@ void AtomVecEllipsoid::grow(int n)
   v = memory->grow(atom->v,nmax,3,"atom:v");
   f = memory->grow(atom->f,nmax,3,"atom:f");
 
+  shape = memory->grow(atom->shape,nmax,3,"atom:shape");
+  rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
   quat = memory->grow(atom->quat,nmax,4,"atom:quat");
   angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
   torque = memory->grow(atom->torque,nmax,3,"atom:torque");
@@ -92,6 +96,7 @@ void AtomVecEllipsoid::grow_reset()
   tag = atom->tag; type = atom->type;
   mask = atom->mask; image = atom->image;
   x = atom->x; v = atom->v; f = atom->f;
+  shape = atom->shape; rmass = atom->rmass;
   quat = atom->quat; angmom = atom->angmom; torque = atom->torque;
 }
 
@@ -110,6 +115,10 @@ void AtomVecEllipsoid::copy(int i, int j)
   v[j][1] = v[i][1];
   v[j][2] = v[i][2];
 
+  shape[j][0] = shape[i][0];
+  shape[j][1] = shape[i][1];
+  shape[j][2] = shape[i][2];
+  rmass[j] = rmass[i];
   quat[j][0] = quat[i][0];
   quat[j][1] = quat[i][1];
   quat[j][2] = quat[i][2];
@@ -392,6 +401,9 @@ int AtomVecEllipsoid::pack_border(int n, int *list, double *buf,
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
+      buf[m++] = shape[j][0];
+      buf[m++] = shape[j][1];
+      buf[m++] = shape[j][2];
       buf[m++] = quat[j][0];
       buf[m++] = quat[j][1];
       buf[m++] = quat[j][2];
@@ -415,6 +427,9 @@ int AtomVecEllipsoid::pack_border(int n, int *list, double *buf,
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
+      buf[m++] = shape[j][0];
+      buf[m++] = shape[j][1];
+      buf[m++] = shape[j][2];
       buf[m++] = quat[j][0];
       buf[m++] = quat[j][1];
       buf[m++] = quat[j][2];
@@ -442,6 +457,9 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
+      buf[m++] = shape[j][0];
+      buf[m++] = shape[j][1];
+      buf[m++] = shape[j][2];
       buf[m++] = quat[j][0];
       buf[m++] = quat[j][1];
       buf[m++] = quat[j][2];
@@ -472,6 +490,9 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
 	buf[m++] = tag[j];
 	buf[m++] = type[j];
 	buf[m++] = mask[j];
+	buf[m++] = shape[j][0];
+	buf[m++] = shape[j][1];
+	buf[m++] = shape[j][2];
 	buf[m++] = quat[j][0];
 	buf[m++] = quat[j][1];
 	buf[m++] = quat[j][2];
@@ -495,6 +516,9 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
 	buf[m++] = tag[j];
 	buf[m++] = type[j];
 	buf[m++] = mask[j];
+	buf[m++] = shape[j][0];
+	buf[m++] = shape[j][1];
+	buf[m++] = shape[j][2];
 	buf[m++] = quat[j][0];
 	buf[m++] = quat[j][1];
 	buf[m++] = quat[j][2];
@@ -521,11 +545,14 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
 
 int AtomVecEllipsoid::pack_border_one(int i, double *buf)
 {
-  buf[0] = quat[i][0];
+  buf[0] = shape[i][0];
-  buf[1] = quat[i][1];
+  buf[1] = shape[i][1];
-  buf[2] = quat[i][2];
+  buf[2] = shape[i][2];
-  buf[3] = quat[i][3];
+  buf[3] = quat[i][0];
-  return 4;
+  buf[4] = quat[i][1];
+  buf[5] = quat[i][2];
+  buf[6] = quat[i][3];
+  return 7;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -544,6 +571,9 @@ void AtomVecEllipsoid::unpack_border(int n, int first, double *buf)
     tag[i] = static_cast<int> (buf[m++]);
     type[i] = static_cast<int> (buf[m++]);
     mask[i] = static_cast<int> (buf[m++]);
+    shape[i][0] = buf[m++];
+    shape[i][1] = buf[m++];
+    shape[i][2] = buf[m++];
     quat[i][0] = buf[m++];
     quat[i][1] = buf[m++];
     quat[i][2] = buf[m++];
@@ -567,6 +597,9 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf)
     tag[i] = static_cast<int> (buf[m++]);
     type[i] = static_cast<int> (buf[m++]);
     mask[i] = static_cast<int> (buf[m++]);
+    shape[i][0] = buf[m++];
+    shape[i][1] = buf[m++];
+    shape[i][2] = buf[m++];
     quat[i][0] = buf[m++];
     quat[i][1] = buf[m++];
     quat[i][2] = buf[m++];
@@ -584,11 +617,14 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf)
 
 int AtomVecEllipsoid::unpack_border_one(int i, double *buf)
 {
-  quat[i][0] = buf[0];
+  shape[i][0] = buf[0];
-  quat[i][1] = buf[1];
+  shape[i][1] = buf[1];
-  quat[i][2] = buf[2];
+  shape[i][2] = buf[2];
-  quat[i][3] = buf[3];
+  quat[i][0] = buf[3];
-  return 4;
+  quat[i][1] = buf[4];
+  quat[i][2] = buf[5];
+  quat[i][3] = buf[6];
+  return 7;
 }
 
 /* ----------------------------------------------------------------------
@@ -610,6 +646,10 @@ int AtomVecEllipsoid::pack_exchange(int i, double *buf)
   buf[m++] = mask[i];
   buf[m++] = image[i];
 
+  buf[m++] = shape[i][0];
+  buf[m++] = shape[i][1];
+  buf[m++] = shape[i][2];
+  buf[m++] = rmass[i];
   buf[m++] = quat[i][0];
   buf[m++] = quat[i][1];
   buf[m++] = quat[i][2];
@@ -645,6 +685,10 @@ int AtomVecEllipsoid::unpack_exchange(double *buf)
   mask[nlocal] = static_cast<int> (buf[m++]);
   image[nlocal] = static_cast<int> (buf[m++]);
 
+  shape[nlocal][0] = buf[m++];
+  shape[nlocal][1] = buf[m++];
+  shape[nlocal][2] = buf[m++];
+  rmass[nlocal] = buf[m++];
   quat[nlocal][0] = buf[m++];
   quat[nlocal][1] = buf[m++];
   quat[nlocal][2] = buf[m++];
@@ -672,7 +716,7 @@ int AtomVecEllipsoid::size_restart()
   int i;
 
   int nlocal = atom->nlocal;
-  int n = 18 * nlocal;
+  int n = 22 * nlocal;
 
   if (atom->nextra_restart)
     for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@@ -702,6 +746,10 @@ int AtomVecEllipsoid::pack_restart(int i, double *buf)
   buf[m++] = v[i][1];
   buf[m++] = v[i][2];
 
+  buf[m++] = shape[i][0];
+  buf[m++] = shape[i][1];
+  buf[m++] = shape[i][2];
+  buf[m++] = rmass[i];
   buf[m++] = quat[i][0];
   buf[m++] = quat[i][1];
   buf[m++] = quat[i][2];
@@ -743,6 +791,10 @@ int AtomVecEllipsoid::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
 
+  shape[nlocal][0] = buf[m++];
+  shape[nlocal][1] = buf[m++];
+  shape[nlocal][2] = buf[m++];
+  rmass[nlocal] = buf[m++];
   quat[nlocal][0] = buf[m++];
   quat[nlocal][1] = buf[m++];
   quat[nlocal][2] = buf[m++];
@@ -781,6 +833,12 @@ void AtomVecEllipsoid::create_atom(int itype, double *coord)
   v[nlocal][0] = 0.0;
   v[nlocal][1] = 0.0;
   v[nlocal][2] = 0.0;
+  
+  shape[nlocal][0] = 0.5;
+  shape[nlocal][1] = 0.5;
+  shape[nlocal][2] = 0.5;
+  rmass[nlocal] = 4.0*PI/3.0 *
+    shape[nlocal][0]*shape[nlocal][1]*shape[nlocal][2];
 
   quat[nlocal][0] = 1.0;
   quat[nlocal][1] = 0.0;
@@ -811,14 +869,36 @@ void AtomVecEllipsoid::data_atom(double *coord, int imagetmp, char **values)
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one("Invalid atom type in Atoms section of data file");
 
+  shape[nlocal][0] = 0.5 * atof(values[2]);
+  shape[nlocal][1] = 0.5 * atof(values[3]);
+  shape[nlocal][2] = 0.5 * atof(values[4]);
+  if (shape[nlocal][0] < 0.0 || shape[nlocal][1] < 0.0 || 
+      shape[nlocal][2] < 0.0)
+    error->one("Invalid shape in Atoms section of data file");
+  if (shape[nlocal][0] > 0.0 || shape[nlocal][1] > 0.0 || 
+      shape[nlocal][2] > 0.0) {
+    if (shape[nlocal][0] == 0.0 || shape[nlocal][1] == 0.0 || 
+	shape[nlocal][2] == 0.0)
+      error->one("Invalid shape in Atoms section of data file");
+  }
+
+  double density = atof(values[5]);
+  if (density <= 0.0)
+    error->one("Invalid density in Atoms section of data file");
+
+  if (shape[nlocal][0] == 0.0) rmass[nlocal] = density;
+  else 
+    rmass[nlocal] = 4.0*PI/3.0 *
+      shape[nlocal][0]*shape[nlocal][1]*shape[nlocal][2] * density;
+
   x[nlocal][0] = coord[0];
   x[nlocal][1] = coord[1];
   x[nlocal][2] = coord[2];
 
-  quat[nlocal][0] = atof(values[5]);
+  quat[nlocal][0] = atof(values[9]);
-  quat[nlocal][1] = atof(values[6]);
+  quat[nlocal][1] = atof(values[10]);
-  quat[nlocal][2] = atof(values[7]);
+  quat[nlocal][2] = atof(values[11]);
-  quat[nlocal][3] = atof(values[8]);
+  quat[nlocal][3] = atof(values[12]);
   MathExtra::normalize4(quat[nlocal]);
 
   image[nlocal] = imagetmp;
@@ -841,13 +921,35 @@ void AtomVecEllipsoid::data_atom(double *coord, int imagetmp, char **values)
 
 int AtomVecEllipsoid::data_atom_hybrid(int nlocal, char **values)
 {
-  quat[nlocal][0] = atof(values[0]);
+  shape[nlocal][0] = 0.5 * atof(values[0]);
-  quat[nlocal][1] = atof(values[1]);
+  shape[nlocal][1] = 0.5 * atof(values[1]);
-  quat[nlocal][2] = atof(values[2]);
+  shape[nlocal][2] = 0.5 * atof(values[2]);
-  quat[nlocal][3] = atof(values[3]);
+  if (shape[nlocal][0] < 0.0 || shape[nlocal][1] < 0.0 || 
+      shape[nlocal][2] < 0.0)
+    error->one("Invalid shape in Atoms section of data file");
+  if (shape[nlocal][0] > 0.0 || shape[nlocal][1] > 0.0 || 
+      shape[nlocal][2] > 0.0) {
+    if (shape[nlocal][0] == 0.0 || shape[nlocal][1] == 0.0 || 
+	shape[nlocal][2] == 0.0)
+      error->one("Invalid shape in Atoms section of data file");
+  }
+
+  double density = atof(values[3]);
+  if (density <= 0.0)
+    error->one("Invalid density in Atoms section of data file");
+
+  if (shape[nlocal][0] == 0.0) rmass[nlocal] = density;
+  else 
+    rmass[nlocal] = 4.0*PI/3.0 *
+      shape[nlocal][0]*shape[nlocal][1]*shape[nlocal][2] * density;
+
+  quat[nlocal][0] = atof(values[4]);
+  quat[nlocal][1] = atof(values[5]);
+  quat[nlocal][2] = atof(values[6]);
+  quat[nlocal][3] = atof(values[7]);
   MathExtra::normalize4(quat[nlocal]);
 
-  return 4;
+  return 8;
 }
 
 /* ----------------------------------------------------------------------
@@ -892,6 +994,8 @@ bigint AtomVecEllipsoid::memory_usage()
   if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
   if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
 
+  if (atom->memcheck("shape")) bytes += memory->usage(shape,nmax,3);
+  if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
   if (atom->memcheck("quat")) bytes += memory->usage(quat,nmax,4);
   if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
   if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);

src/ASPHERE/atom_vec_ellipsoid.h

@@ -60,8 +60,10 @@ class AtomVecEllipsoid : public AtomVec {
   bigint memory_usage();
 
  private:
+  double PI;
   int *tag,*type,*mask,*image;
   double **x,**v,**f;
+  double **shape,*density,*rmass;
   double **angmom,**torque,**quat;
 };
 

src/ASPHERE/compute_erotate_asphere.cpp

@@ -34,24 +34,10 @@ ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   extscalar = 1;
 
-  memory->create(inertia,atom->ntypes+1,3,"fix_temp_sphere:inertia");
+  // error check
 
-  // error checks
+  if (!atom->ellipsoid_flag)
-
+    error->all("Compute erotate/asphere requires atom style ellipsoid");
-  if (!atom->angmom_flag || !atom->quat_flag ||
-      !atom->avec->shape_type)
-    error->all("Compute erotate/asphere requires atom attributes "
-	       "angmom, quat, shape");
-  if (atom->radius_flag || atom->rmass_flag)
-    error->all("Compute erotate/asphere cannot be used with atom attributes "
-	       "diameter or rmass");
-}
-
-/* ---------------------------------------------------------------------- */
-
-ComputeERotateAsphere::~ComputeERotateAsphere()
-{
-  memory->destroy(inertia);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -62,77 +48,62 @@ void ComputeERotateAsphere::init()
   // no point particles allowed, spherical is OK
 
   double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit)
-      if (shape[type[i]][0] == 0.0)
+      if (shape[i][0] == 0.0)
 	error->one("Compute erotate/asphere requires extended particles");
 
   pfactor = 0.5 * force->mvv2e;
-  calculate_inertia();
 }
 
 /* ---------------------------------------------------------------------- */
 
 double ComputeERotateAsphere::compute_scalar()
 {
-  int i,itype;
-
   invoked_scalar = update->ntimestep;
 
   double **quat = atom->quat;
   double **angmom = atom->angmom;
+  double **shape = atom->shape;
+  double *rmass = atom->rmass;
   int *mask = atom->mask;
-  int *type = atom->type;
   int nlocal = atom->nlocal;
 
   // sum rotational energy for each particle
   // no point particles since divide by inertia
 
-  double wbody[3];
+  double wbody[3],inertia[3];
   double rot[3][3];
   double erotate = 0.0;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      itype = type[i];
+
+      // principal moments of inertia
+
+      inertia[0] = rmass[i] * 
+	(shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[1] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[2] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
 
       // wbody = angular velocity in body frame
 
       MathExtra::quat_to_mat(quat[i],rot);
       MathExtra::transpose_times_column3(rot,angmom[i],wbody);
-      wbody[0] /= inertia[itype][0];
+      wbody[0] /= inertia[0];
-      wbody[1] /= inertia[itype][1];
+      wbody[1] /= inertia[1];
-      wbody[2] /= inertia[itype][2];
+      wbody[2] /= inertia[2];
       
-      erotate += inertia[itype][0]*wbody[0]*wbody[0]+
+      erotate += inertia[0]*wbody[0]*wbody[0] +
-	inertia[itype][1]*wbody[1]*wbody[1]+
+	inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
-	inertia[itype][2]*wbody[2]*wbody[2];
     }
 
   MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   scalar *= pfactor;
   return scalar;
 }
-
-/* ----------------------------------------------------------------------
-   principal moments of inertia for ellipsoids
-------------------------------------------------------------------------- */
-
-void ComputeERotateAsphere::calculate_inertia()
-{
-  double *mass = atom->mass;
-  double **shape = atom->shape;
-
-  for (int i = 1; i <= atom->ntypes; i++) {
-    inertia[i][0] = mass[i] * 
-      (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
-    inertia[i][1] = mass[i] * 
-      (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
-    inertia[i][2] = mass[i] * 
-      (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
-  }
-}

src/ASPHERE/compute_erotate_asphere.h

@@ -27,15 +27,11 @@ namespace LAMMPS_NS {
 class ComputeERotateAsphere : public Compute {
  public:
   ComputeERotateAsphere(class LAMMPS *, int, char **);
-  ~ComputeERotateAsphere();
   void init();
   double compute_scalar();
 
  private:
   double pfactor;
-  double **inertia;
-
-  void calculate_inertia();
 };
 
 }

src/ASPHERE/compute_temp_asphere.cpp

@@ -55,17 +55,11 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
   }
 
   vector = new double[6];
-  memory->create(inertia,atom->ntypes+1,3,"fix_temp_sphere:inertia");
 
   // error checks
 
-  if (!atom->angmom_flag || !atom->quat_flag ||
+  if (!atom->ellipsoid_flag) 
-      !atom->avec->shape_type)
+    error->all("Compute temp/asphere requires atom style ellipsoid");
-    error->all("Compute temp/asphere requires atom attributes "
-	       "angmom, quat, shape");
-  if (atom->radius_flag || atom->rmass_flag)
-    error->all("Compute temp/asphere cannot be used with atom attributes "
-	       "diameter or rmass");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -74,7 +68,6 @@ ComputeTempAsphere::~ComputeTempAsphere()
 {
   delete [] id_bias;
   delete [] vector;
-  memory->destroy(inertia);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -85,13 +78,12 @@ void ComputeTempAsphere::init()
   // no point particles allowed, spherical is OK
 
   double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit)
-      if (shape[type[i]][0] == 0.0)
+      if (shape[i][0] == 0.0)
 	error->one("Compute temp/asphere requires extended particles");
 
   if (tempbias) {
@@ -113,8 +105,6 @@ void ComputeTempAsphere::init()
   for (int i = 0; i < modify->nfix; i++)
     fix_dof += modify->fix[i]->dof(igroup);
   dof_compute();
-
-  calculate_inertia();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -164,13 +154,12 @@ double ComputeTempAsphere::compute_scalar()
   double **v = atom->v;
   double **quat = atom->quat;
   double **angmom = atom->angmom;
-  double *mass = atom->mass;
+  double **shape = atom->shape;
-  int *type = atom->type;
+  double *rmass = atom->rmass;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  int itype;
+  double wbody[3],inertia[3];
-  double wbody[3];
   double rot[3][3];
   double t = 0.0;
   
@@ -180,20 +169,27 @@ double ComputeTempAsphere::compute_scalar()
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
 
-      itype = type[i];
+      t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
-      t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[itype];
+
+      // principal moments of inertia
+
+      inertia[0] = rmass[i] * 
+	(shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[1] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[2] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
 
       // wbody = angular velocity in body frame
       
       MathExtra::quat_to_mat(quat[i],rot);
       MathExtra::transpose_times_column3(rot,angmom[i],wbody);
-      wbody[0] /= inertia[itype][0];
+      wbody[0] /= inertia[0];
-      wbody[1] /= inertia[itype][1];
+      wbody[1] /= inertia[1];
-      wbody[2] /= inertia[itype][2];
+      wbody[2] /= inertia[2];
 
-      t += inertia[itype][0]*wbody[0]*wbody[0] +
+      t += inertia[0]*wbody[0]*wbody[0] +
-	inertia[itype][1]*wbody[1]*wbody[1] +
+	inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
-	inertia[itype][2]*wbody[2]*wbody[2];
     }
 
   if (tempbias) tbias->restore_bias_all();
@@ -220,13 +216,12 @@ void ComputeTempAsphere::compute_vector()
   double **v = atom->v;
   double **quat = atom->quat;
   double **angmom = atom->angmom;
-  double *mass = atom->mass;
+  double **shape = atom->shape;
-  int *type = atom->type;
+  double *rmass = atom->rmass;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  int itype;
+  double wbody[3],inertia[3];
-  double wbody[3];
   double rot[3][3];
   double massone,t[6];
   for (i = 0; i < 6; i++) t[i] = 0.0;
@@ -236,8 +231,7 @@ void ComputeTempAsphere::compute_vector()
 
       // translational kinetic energy
 
-      itype = type[i];
+      massone = rmass[i];
-      massone = mass[itype];
       t[0] += massone * v[i][0]*v[i][0];
       t[1] += massone * v[i][1]*v[i][1];
       t[2] += massone * v[i][2]*v[i][2];
@@ -245,22 +239,31 @@ void ComputeTempAsphere::compute_vector()
       t[4] += massone * v[i][0]*v[i][2];
       t[5] += massone * v[i][1]*v[i][2];
 
+      // principal moments of inertia
+
+      inertia[0] = rmass[i] * 
+	(shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[1] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[2] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+
       // wbody = angular velocity in body frame
 
       MathExtra::quat_to_mat(quat[i],rot);
       MathExtra::transpose_times_column3(rot,angmom[i],wbody);
-      wbody[0] /= inertia[itype][0];
+      wbody[0] /= inertia[0];
-      wbody[1] /= inertia[itype][1];
+      wbody[1] /= inertia[1];
-      wbody[2] /= inertia[itype][2];
+      wbody[2] /= inertia[2];
 
       // rotational kinetic energy
 
-      t[0] += inertia[itype][0]*wbody[0]*wbody[0];
+      t[0] += inertia[0]*wbody[0]*wbody[0];
-      t[1] += inertia[itype][1]*wbody[1]*wbody[1];
+      t[1] += inertia[1]*wbody[1]*wbody[1];
-      t[2] += inertia[itype][2]*wbody[2]*wbody[2];
+      t[2] += inertia[2]*wbody[2]*wbody[2];
-      t[3] += inertia[itype][0]*wbody[0]*wbody[1];
+      t[3] += inertia[0]*wbody[0]*wbody[1];
-      t[4] += inertia[itype][1]*wbody[0]*wbody[2];
+      t[4] += inertia[1]*wbody[0]*wbody[2];
-      t[5] += inertia[itype][2]*wbody[1]*wbody[2];
+      t[5] += inertia[2]*wbody[1]*wbody[2];
     }
 
   if (tempbias) tbias->restore_bias_all();
@@ -269,25 +272,6 @@ void ComputeTempAsphere::compute_vector()
   for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
 }
 
-/* ----------------------------------------------------------------------
-   principal moments of inertia for ellipsoids
-------------------------------------------------------------------------- */
-
-void ComputeTempAsphere::calculate_inertia()
-{
-  double *mass = atom->mass;
-  double **shape = atom->shape;
-
-  for (int i = 1; i <= atom->ntypes; i++) {
-    inertia[i][0] = mass[i] * 
-      (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
-    inertia[i][1] = mass[i] * 
-      (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
-    inertia[i][2] = mass[i] * 
-      (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
-  }
-}
-
 /* ----------------------------------------------------------------------
    remove velocity bias from atom I to leave thermal velocity
 ------------------------------------------------------------------------- */

src/ASPHERE/compute_temp_asphere.h

@@ -38,12 +38,10 @@ class ComputeTempAsphere : public Compute {
  private:
   int fix_dof;
   double tfactor;
-  double **inertia;
   char *id_bias;
   Compute *tbias;     // ptr to additional bias compute
 
   void dof_compute();
-  void calculate_inertia();
 };
 
 }

src/ASPHERE/fix_nh_asphere.cpp

@@ -33,22 +33,8 @@ using namespace LAMMPS_NS;
 FixNHAsphere::FixNHAsphere(LAMMPS *lmp, int narg, char **arg) :
   FixNH(lmp, narg, arg)
 {
-  memory->create(inertia,atom->ntypes+1,3,"fix_nvt_asphere:inertia");
+  if (!atom->ellipsoid_flag)
-
+    error->all("Fix nvt/nph/npt asphere requires atom style ellipsoid");
-  if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag ||
-      !atom->avec->shape_type)
-    error->all("Fix nvt/nph/npt asphere requires atom attributes "
-	       "quat, angmom, torque, shape");
-  if (atom->radius_flag || atom->rmass_flag)
-    error->all("Fix nvt/nph/npt asphere cannot be used with atom attributes "
-	       "diameter or rmass");
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixNHAsphere::~FixNHAsphere()
-{
-  memory->destroy(inertia);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -59,17 +45,15 @@ void FixNHAsphere::init()
   // no point particles allowed, spherical is OK
 
   double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit)
-      if (shape[type[i]][0] == 0.0)
+      if (shape[i][0] == 0.0)
 	error->one("Fix nvt/nph/npt asphere requires extended particles");
 
   FixNH::init();
-  calculate_inertia();
 }
 
 /* ----------------------------------------------------------------------
@@ -149,25 +133,6 @@ void FixNHAsphere::omega_from_mq(double *q, double *m, double *inertia,
   MathExtra::times_column3(rot,wbody,w);
 }
 
-/* ----------------------------------------------------------------------
-   principal moments of inertia for ellipsoids
-------------------------------------------------------------------------- */
-
-void FixNHAsphere::calculate_inertia()
-{
-  double *mass = atom->mass;
-  double **shape = atom->shape;
-  
-  for (int i = 1; i <= atom->ntypes; i++) {
-    inertia[i][0] = 0.2*mass[i] *
-      (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]);
-    inertia[i][1] = 0.2*mass[i] *
-      (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]);
-    inertia[i][2] = 0.2*mass[i] * 
-      (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]);
-  }
-}
-
 /* ----------------------------------------------------------------------
    perform half-step update of angular momentum
 -----------------------------------------------------------------------*/
@@ -207,7 +172,8 @@ void FixNHAsphere::nve_x()
 
   double **quat = atom->quat;
   double **angmom = atom->angmom;
-  int *type = atom->type;
+  double **shape = atom->shape;
+  double *rmass = atom->rmass;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
@@ -218,11 +184,21 @@ void FixNHAsphere::nve_x()
 
   // update quaternion a full step via Richardson iteration
   // returns new normalized quaternion
+  // principal moments of inertia
 
-  for (int i = 0; i < nlocal; i++) {    
+  double inertia[3];
-    if (mask[i] & groupbit)
+
-      richardson(quat[i],angmom[i],inertia[type[i]]);
+  for (int i = 0; i < nlocal; i++)
-  }
+    if (mask[i] & groupbit) {
+      inertia[0] = rmass[i] * 
+	(shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[1] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[2] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+      
+      richardson(quat[i],angmom[i],inertia);
+    }
 }
 
 /* ----------------------------------------------------------------------

src/ASPHERE/fix_nh_asphere.h

@@ -21,7 +21,7 @@ namespace LAMMPS_NS {
 class FixNHAsphere : public FixNH {
  public:
   FixNHAsphere(class LAMMPS *, int, char **);
-  virtual ~FixNHAsphere();
+  virtual ~FixNHAsphere() {}
   void init();
 
  protected:

src/ASPHERE/fix_nve_asphere.cpp

@@ -34,24 +34,10 @@ using namespace LAMMPS_NS;
 FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) : 
   FixNVE(lmp, narg, arg)
 {
-  memory->create(inertia,atom->ntypes+1,3,"fix_nve_asphere:inertia");
-
   // error checks
 
-  if (!atom->angmom_flag || !atom->quat_flag || !atom->torque_flag ||
+  if (!atom->ellipsoid_flag)
-      !atom->avec->shape_type)
+    error->all("Fix nve/asphere requires atom style ellipsoid");
-    error->all("Fix nve/asphere requires atom attributes "
-	       "angmom, quat, torque, shape");
-  if (atom->radius_flag || atom->rmass_flag)
-    error->all("Fix nve/asphere cannot be used with atom attributes "
-	       "diameter or rmass");
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixNVEAsphere::~FixNVEAsphere()
-{
-  memory->destroy(inertia);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -62,17 +48,15 @@ void FixNVEAsphere::init()
   // no point particles allowed, spherical is OK
 
   double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit)
-      if (shape[type[i]][0] == 0.0)
+      if (shape[i][0] == 0.0)
 	error->one("Fix nve/asphere requires extended particles");
 
   FixNVE::init();
-  calculate_inertia();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -80,6 +64,7 @@ void FixNVEAsphere::init()
 void FixNVEAsphere::initial_integrate(int vflag)
 {
   double dtfm;
+  double inertia[3];
 
   double **x = atom->x;
   double **v = atom->v;
@@ -87,8 +72,8 @@ void FixNVEAsphere::initial_integrate(int vflag)
   double **quat = atom->quat;
   double **angmom = atom->angmom;
   double **torque = atom->torque;
-  double *mass = atom->mass;
+  double **shape = atom->shape;
-  int *type = atom->type;
+  double *rmass = atom->rmass;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
@@ -99,7 +84,7 @@ void FixNVEAsphere::initial_integrate(int vflag)
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      dtfm = dtf / mass[type[i]];
+      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];
@@ -115,7 +100,16 @@ void FixNVEAsphere::initial_integrate(int vflag)
       angmom[i][1] += dtf * torque[i][1];
       angmom[i][2] += dtf * torque[i][2];
 
-      richardson(quat[i],angmom[i],inertia[type[i]]);
+      // principal moments of inertia
+
+      inertia[0] = rmass[i] * 
+	(shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[1] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
+      inertia[2] = rmass[i] * 
+	(shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+
+      richardson(quat[i],angmom[i],inertia);
     }
 }
 
@@ -129,15 +123,14 @@ void FixNVEAsphere::final_integrate()
   double **f = atom->f;
   double **angmom = atom->angmom;
   double **torque = atom->torque;
-  double *mass = atom->mass;
+  double *rmass = atom->rmass;
-  int *type = atom->type;
   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 / mass[type[i]];
+      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];
@@ -224,22 +217,3 @@ void FixNVEAsphere::omega_from_mq(double *q, double *m, double *moments,
   wbody[2] /= moments[2];
   MathExtra::times_column3(rot,wbody,w);
 }
-
-/* ----------------------------------------------------------------------
-   principal moments of inertia for ellipsoids
-------------------------------------------------------------------------- */
-
-void FixNVEAsphere::calculate_inertia()
-{
-  double *mass = atom->mass;
-  double **shape = atom->shape;
-  
-  for (int i = 1; i <= atom->ntypes; i++) {
-    inertia[i][0] = 0.2*mass[i] *
-      (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]);
-    inertia[i][1] = 0.2*mass[i] *
-      (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]);
-    inertia[i][2] = 0.2*mass[i] * 
-      (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]);
-  }
-}

src/ASPHERE/fix_nve_asphere.h

@@ -27,18 +27,15 @@ namespace LAMMPS_NS {
 class FixNVEAsphere : public FixNVE {
  public:
   FixNVEAsphere(class LAMMPS *, int, char **);
-  ~FixNVEAsphere();
   void init();
   void initial_integrate(int);
   void final_integrate();
 
  private:
   double dtq;
-  double **inertia;
 
   void richardson(double *, double *, double *);
   void omega_from_mq(double *, double *, double *, double *);
-  void calculate_inertia();
 };
 
 }

src/ASPHERE/pair_gayberne.cpp

@@ -58,7 +58,8 @@ PairGayBerne::~PairGayBerne()
     memory->destroy(form);
     memory->destroy(epsilon);
     memory->destroy(sigma);
-    memory->destroy(shape);
+    memory->destroy(shape1);
+    memory->destroy(shape2);
     memory->destroy(well);
     memory->destroy(cut);
     memory->destroy(lj1);
@@ -109,7 +110,7 @@ void PairGayBerne::compute(int eflag, int vflag)
       MathExtra::quat_to_mat_trans(quat[i],a1);
       MathExtra::diag_times3(well[itype],a1,temp);
       MathExtra::transpose_times3(a1,temp,b1);
-      MathExtra::diag_times3(shape[itype],a1,temp);
+      MathExtra::diag_times3(shape2[itype],a1,temp);
       MathExtra::transpose_times3(a1,temp,g1);
     }
 
@@ -153,7 +154,7 @@ void PairGayBerne::compute(int eflag, int vflag)
 	  MathExtra::quat_to_mat_trans(quat[j],a2);
 	  MathExtra::diag_times3(well[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,b2);
-	  MathExtra::diag_times3(shape[jtype],a2,temp);
+	  MathExtra::diag_times3(shape2[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,g2);
 	  one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor);
 	  ttor[0] = ttor[1] = ttor[2] = 0.0;
@@ -168,7 +169,7 @@ void PairGayBerne::compute(int eflag, int vflag)
 	  MathExtra::quat_to_mat_trans(quat[j],a2);
 	  MathExtra::diag_times3(well[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,b2);
-	  MathExtra::diag_times3(shape[jtype],a2,temp);
+	  MathExtra::diag_times3(shape2[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,g2);
 	  one_eng = gayberne_analytic(i,j,a1,a2,b1,b2,g1,g2,r12,rsq,
 				      fforce,ttor,rtor);
@@ -232,7 +233,8 @@ void PairGayBerne::allocate()
   memory->create(form,n+1,n+1,"pair:form");
   memory->create(epsilon,n+1,n+1,"pair:epsilon");
   memory->create(sigma,n+1,n+1,"pair:sigma");
-  memory->create(shape,n+1,3,"pair:shape");
+  memory->create(shape1,n+1,3,"pair:shape1");
+  memory->create(shape2,n+1,3,"pair:shape2");
   memory->create(well,n+1,3,"pair:well");
   memory->create(cut,n+1,n+1,"pair:cut");
   memory->create(lj1,n+1,n+1,"pair:lj1");
@@ -330,22 +332,23 @@ void PairGayBerne::coeff(int narg, char **arg)
 
 void PairGayBerne::init_style()
 {
-  if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type)
+  if (!atom->ellipsoid_flag)
-    error->all("Pair gayberne requires atom attributes quat, torque, shape");
+    error->all("Pair gayberne requires atom style ellipsoid");
-  if (atom->radius_flag)
-    error->all("Pair gayberne cannot be used with atom attribute diameter");
 
   int irequest = neighbor->request(this);
 
   // per-type shape precalculations
+  // require that atom shapes are identical within each type
 
   for (int i = 1; i <= atom->ntypes; i++) {
+    if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2]))
+      error->all("Pair gayberne requires atoms with same type have same shape");
     if (setwell[i]) {
-      double *one = atom->shape[i];
+      shape2[i][0] = shape1[i][0]*shape1[i][0];
-      shape[i][0] = one[0]*one[0];
+      shape2[i][1] = shape1[i][1]*shape1[i][1];
-      shape[i][1] = one[1]*one[1];
+      shape2[i][2] = shape1[i][2]*shape1[i][2];
-      shape[i][2] = one[2]*one[2];
+      lshape[i] = (shape1[i][0]*shape1[i][1]+shape1[i][2]*shape1[i][2]) * 
-      lshape[i] = (one[0]*one[1]+one[2]*one[2])*sqrt(one[0]*one[1]);
+	sqrt(shape1[i][0]*shape1[i][1]);
     }
   }
 }
@@ -377,14 +380,14 @@ double PairGayBerne::init_one(int i, int j)
   } else offset[i][j] = 0.0;
 
   int ishape = 0;
-  if (atom->shape[i][0] != atom->shape[i][1] || 
+  if (shape1[i][0] != shape1[i][1] || 
-      atom->shape[i][0] != atom->shape[i][2] || 
+      shape1[i][0] != shape1[i][2] || 
-      atom->shape[i][1] != atom->shape[i][2]) ishape = 1;
+      shape1[i][1] != shape1[i][2]) ishape = 1;
   if (setwell[i] == 1) ishape = 1;
   int jshape = 0;
-  if (atom->shape[j][0] != atom->shape[j][1] || 
+  if (shape1[j][0] != shape1[j][1] || 
-      atom->shape[j][0] != atom->shape[j][2] || 
+      shape1[j][0] != shape1[j][2] || 
-      atom->shape[j][1] != atom->shape[j][2]) jshape = 1;
+      shape1[j][1] != shape1[j][2]) jshape = 1;
   if (setwell[j] == 1) jshape = 1;
   
   if (ishape == 0 && jshape == 0)
@@ -640,7 +643,7 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
 
   double deta[3];
   deta[0] = deta[1] = deta[2] = 0.0;
-  compute_eta_torque(g12,a1,shape[type[i]],temp);
+  compute_eta_torque(g12,a1,shape2[type[i]],temp);
   temp1 = -eta*upsilon;
   for (int m = 0; m < 3; m++) {
     for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y];
@@ -655,7 +658,7 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
   double deta2[3];
   if (newton_pair || j < nlocal) {
     deta2[0] = deta2[1] = deta2[2] = 0.0;
-    compute_eta_torque(g12,a2,shape[type[j]],temp);
+    compute_eta_torque(g12,a2,shape2[type[j]],temp);
     for (int m = 0; m < 3; m++) {
       for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y];
       MathExtra::cross3(a2[m],tempv,tempv2);
@@ -707,9 +710,9 @@ double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3],
   // compute distance of closest approach
 
   double g12[3][3];
-  g12[0][0] = g1[0][0]+shape[type[j]][0];  
+  g12[0][0] = g1[0][0]+shape2[type[j]][0];  
-  g12[1][1] = g1[1][1]+shape[type[j]][0];  
+  g12[1][1] = g1[1][1]+shape2[type[j]][0];  
-  g12[2][2] = g1[2][2]+shape[type[j]][0];  
+  g12[2][2] = g1[2][2]+shape2[type[j]][0];  
   g12[0][1] = g1[0][1]; g12[1][0] = g1[1][0];
   g12[0][2] = g1[0][2]; g12[2][0] = g1[2][0];
   g12[1][2] = g1[1][2]; g12[2][1] = g1[2][1];
@@ -809,7 +812,7 @@ double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3],
 
   double deta[3];
   deta[0] = deta[1] = deta[2] = 0.0;
-  compute_eta_torque(g12,a1,shape[type[i]],temp);
+  compute_eta_torque(g12,a1,shape2[type[i]],temp);
   temp1 = -eta*upsilon;
   for (int m = 0; m < 3; m++) {
     for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y];

src/ASPHERE/pair_gayberne.h

@@ -43,7 +43,8 @@ class PairGayBerne : public Pair {
   double **cut;
 
   double gamma,upsilon,mu;   // Gay-Berne parameters
-  double **shape;            // radii in x, y and z SQUARED
+  double **shape1;           // per-type radii in x, y and z
+  double **shape2;           // per-type radii in x, y and z SQUARED
   double *lshape;            // precalculation based on the shape
   double **well;             // well depth scaling along each axis ^ -1.0/mu
   double **epsilon,**sigma;  // epsilon and sigma values for atom-type pairs

src/ASPHERE/pair_resquared.cpp

@@ -64,6 +64,7 @@ PairRESquared::~PairRESquared()
     memory->destroy(form);
     memory->destroy(epsilon);
     memory->destroy(sigma);
+    memory->destroy(shape1);
     memory->destroy(shape2);
     memory->destroy(well);
     memory->destroy(cut);
@@ -224,6 +225,7 @@ void PairRESquared::allocate()
   memory->create(form,n+1,n+1,"pair:form");
   memory->create(epsilon,n+1,n+1,"pair:epsilon");
   memory->create(sigma,n+1,n+1,"pair:sigma");
+  memory->create(shape1,n+1,3,"pair:shape1");
   memory->create(shape2,n+1,3,"pair:shape2");
   memory->create(well,n+1,3,"pair:well");
   memory->create(cut,n+1,n+1,"pair:cut");
@@ -319,22 +321,22 @@ void PairRESquared::coeff(int narg, char **arg)
 
 void PairRESquared::init_style()
 {
-  if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type)
+  if (!atom->ellipsoid_flag)
-    error->all("Pair resquared requires atom attributes quat, torque, shape");
+    error->all("Pair resquared requires atom style ellipsoid");
-  if (atom->radius_flag)
-    error->all("Pair resquared cannot be used with atom attribute diameter");
 
   int irequest = neighbor->request(this);
 
   // per-type shape precalculations
+  // require that atom shapes are identical within each type
 
   for (int i = 1; i <= atom->ntypes; i++) {
+    if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2]))
+      error->all("Pair gayberne requires atoms with same type have same shape");
     if (setwell[i]) {
-      double *one = atom->shape[i];
+      shape2[i][0] = shape1[i][0]*shape1[i][0];
-      shape2[i][0] = one[0]*one[0];
+      shape2[i][1] = shape1[i][1]*shape1[i][1];
-      shape2[i][1] = one[1]*one[1];
+      shape2[i][2] = shape1[i][2]*shape1[i][2];
-      shape2[i][2] = one[2]*one[2];
+      lshape[i] = shape1[i][0]*shape1[i][1]*shape1[i][2];
-      lshape[i] = one[0]*one[1]*one[2];
     }
   }
 }
@@ -345,16 +347,14 @@ void PairRESquared::init_style()
 
 double PairRESquared::init_one(int i, int j)
 {
-  double **shape = atom->shape;
-  
   if (setwell[i] == 0 || setwell[j] == 0)
     error->all("Pair resquared epsilon a,b,c coeffs are not all set");
 
   int ishape = 0;
-  if (shape[i][0] != 0.0 && shape[i][1] != 0.0 && shape[i][2] != 0.0) 
+  if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0) 
     ishape = 1;
   int jshape = 0;
-  if (shape[j][0] != 0.0 && shape[j][1] != 0.0 && shape[j][2] != 0.0) 
+  if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0) 
     jshape = 1;
   
   if (ishape == 0 && jshape == 0) {
@@ -546,7 +546,6 @@ double PairRESquared::resquared_analytic(const int i, const int j,
                                          double *rtor)
 {
   int *type = atom->type;
-  double **shape = atom->shape;
     
   // pair computations for energy, force, torque
 
@@ -649,16 +648,16 @@ double PairRESquared::resquared_analytic(const int i, const int j,
   tprod = eta*chi*sigh;
 
   double stemp = h12/2.0;
-  Ua = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)*
+  Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
-       (shape[type[i]][2]+stemp)*(shape[type[j]][0]+stemp)*
+       (shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)*
-       (shape[type[j]][1]+stemp)*(shape[type[j]][2]+stemp);
+       (shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp);
   Ua = (1.0+3.0*tprod)*sprod/Ua;
   Ua = epsilon[type[i]][type[j]]*Ua/-36.0;
 
   stemp = h12/cr60;
-  Ur = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)*
+  Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
-       (shape[type[i]][2]+stemp)*(shape[type[j]][0]+stemp)*
+       (shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)*
-       (shape[type[j]][1]+stemp)*(shape[type[j]][2]+stemp);
+       (shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp);
   Ur = (1.0+b_alpha*tprod)*sprod/Ur;
   Ur = epsilon[type[i]][type[j]]*Ur*pow(sigh,6.0)/2025.0;
 
@@ -704,22 +703,22 @@ double PairRESquared::resquared_analytic(const int i, const int j,
   spr[1] = 0.5*sigma12p3*s[1];
   spr[2] = 0.5*sigma12p3*s[2];
 
-  stemp = 1.0/(shape[type[i]][0]*2.0+h12)+
+  stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+
-          1.0/(shape[type[i]][1]*2.0+h12)+
+          1.0/(shape1[type[i]][1]*2.0+h12)+
-          1.0/(shape[type[i]][2]*2.0+h12)+
+          1.0/(shape1[type[i]][2]*2.0+h12)+
-          1.0/(shape[type[j]][0]*2.0+h12)+
+          1.0/(shape1[type[j]][0]*2.0+h12)+
-          1.0/(shape[type[j]][1]*2.0+h12)+
+          1.0/(shape1[type[j]][1]*2.0+h12)+
-          1.0/(shape[type[j]][2]*2.0+h12);
+          1.0/(shape1[type[j]][2]*2.0+h12);
   hsec = h12+3.0*sec;
   dspu = 1.0/h12-1.0/hsec+stemp;
   pbsu = 3.0*sigma[type[i]][type[j]]/hsec;
   
-  stemp = 1.0/(shape[type[i]][0]*cr60+h12)+
+  stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+
-          1.0/(shape[type[i]][1]*cr60+h12)+
+          1.0/(shape1[type[i]][1]*cr60+h12)+
-          1.0/(shape[type[i]][2]*cr60+h12)+
+          1.0/(shape1[type[i]][2]*cr60+h12)+
-          1.0/(shape[type[j]][0]*cr60+h12)+
+          1.0/(shape1[type[j]][0]*cr60+h12)+
-          1.0/(shape[type[j]][1]*cr60+h12)+
+          1.0/(shape1[type[j]][1]*cr60+h12)+
-          1.0/(shape[type[j]][2]*cr60+h12);
+          1.0/(shape1[type[j]][2]*cr60+h12);
   hsec = h12+b_alpha*sec;
   dspr = 7.0/h12-1.0/hsec+stemp;
   pbsr = b_alpha*sigma[type[i]][type[j]]/hsec;
@@ -832,7 +831,6 @@ double PairRESquared::resquared_lj(const int i, const int j,
                                    double *ttor, bool calc_torque)
 {
   int *type = atom->type;
-  double **shape = atom->shape;
     
   // pair computations for energy, force, torque
 
@@ -875,9 +873,9 @@ double PairRESquared::resquared_lj(const int i, const int j,
   double lAtwo[3][3][3];  // A1'*S1^2*wi.lA
   double scorrect[3];
   double half_sigma=sigma[type[i]][type[j]] / 2.0;
-  scorrect[0] = shape[type[i]][0]+half_sigma;
+  scorrect[0] = shape1[type[i]][0]+half_sigma;
-  scorrect[1] = shape[type[i]][1]+half_sigma;
+  scorrect[1] = shape1[type[i]][1]+half_sigma;
-  scorrect[2] = shape[type[i]][2]+half_sigma;
+  scorrect[2] = shape1[type[i]][2]+half_sigma;
   scorrect[0] = scorrect[0] * scorrect[0] / 2.0;
   scorrect[1] = scorrect[1] * scorrect[1] / 2.0;
   scorrect[2] = scorrect[2] * scorrect[2] / 2.0;
@@ -909,14 +907,14 @@ double PairRESquared::resquared_lj(const int i, const int j,
   h12p3 = pow(h12,3.0);
   double sigmap3 = pow(sigma[type[i]][type[j]],3.0);
   double stemp = h12/2.0;
-  Ua = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)*
+  Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
-       (shape[type[i]][2]+stemp)*h12p3/8.0;
+       (shape1[type[i]][2]+stemp)*h12p3/8.0;
   Ua = (1.0+3.0*tprod)*lshape[type[i]]/Ua;
   Ua = epsilon[type[i]][type[j]]*Ua*sigmap3*solv_f_a;
 
   stemp = h12/cr60;
-  Ur = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)*
+  Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
-       (shape[type[i]][2]+stemp)*h12p3/60.0;
+       (shape1[type[i]][2]+stemp)*h12p3/60.0;
   Ur = (1.0+b_alpha*tprod)*lshape[type[i]]/Ur;
   Ur = epsilon[type[i]][type[j]]*Ur*sigmap3*pow(sigh,6.0)*solv_f_r;
 
@@ -931,17 +929,17 @@ double PairRESquared::resquared_lj(const int i, const int j,
   spr[1] = 0.5*sigma12p3*s[1];
   spr[2] = 0.5*sigma12p3*s[2];
 
-  stemp = 1.0/(shape[type[i]][0]*2.0+h12)+
+  stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+
-          1.0/(shape[type[i]][1]*2.0+h12)+
+          1.0/(shape1[type[i]][1]*2.0+h12)+
-          1.0/(shape[type[i]][2]*2.0+h12)+
+          1.0/(shape1[type[i]][2]*2.0+h12)+
           3.0/h12;
   hsec = h12+3.0*sec;
   dspu = 1.0/h12-1.0/hsec+stemp;
   pbsu = 3.0*sigma[type[i]][type[j]]/hsec;
   
-  stemp = 1.0/(shape[type[i]][0]*cr60+h12)+
+  stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+
-          1.0/(shape[type[i]][1]*cr60+h12)+
+          1.0/(shape1[type[i]][1]*cr60+h12)+
-          1.0/(shape[type[i]][2]*cr60+h12)+
+          1.0/(shape1[type[i]][2]*cr60+h12)+
           3.0/h12;
   hsec = h12+b_alpha*sec;
   dspr = 7.0/h12-1.0/hsec+stemp;

src/ASPHERE/pair_resquared.h

@@ -42,7 +42,8 @@ class PairRESquared : public Pair {
   double cut_global;
   double **cut;
 
-  double **shape2;           // radii in x, y and z SQUARED
+  double **shape1;           // per-type radii in x, y and z
+  double **shape2;           // per-type radii in x, y and z SQUARED
   double *lshape;            // product of the radii
   double **well;             // well depth scaling along each axis
   double **epsilon,**sigma;  // epsilon and sigma values for atom-type pairs

src/COLLOID/Install.sh

@@ -2,13 +2,11 @@
 
 if (test $1 = 1) then
 
-  cp atom_vec_colloid.cpp ..
   cp fix_wall_colloid.cpp ..
   cp pair_colloid.cpp ..
   cp pair_lubricate.cpp ..
   cp pair_yukawa_colloid.cpp ..
 
-  cp atom_vec_colloid.h ..
   cp fix_wall_colloid.h ..
   cp pair_colloid.h ..
   cp pair_lubricate.h ..
@@ -16,13 +14,11 @@ if (test $1 = 1) then
 
 elif (test $1 = 0) then
 
-  rm ../atom_vec_colloid.cpp
   rm ../fix_wall_colloid.cpp
   rm ../pair_colloid.cpp
   rm ../pair_lubricate.cpp
   rm ../pair_yukawa_colloid.cpp
 
-  rm ../atom_vec_colloid.h
   rm ../fix_wall_colloid.h
   rm ../pair_colloid.h
   rm ../pair_lubricate.h

src/COLLOID/atom_vec_colloid.cpp

@@ -1,749 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 "lmptype.h"
-#include "stdlib.h"
-#include "atom_vec_colloid.h"
-#include "atom.h"
-#include "force.h"
-#include "domain.h"
-#include "modify.h"
-#include "fix.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-#define DELTA 10000
-
-/* ---------------------------------------------------------------------- */
-
-AtomVecColloid::AtomVecColloid(LAMMPS *lmp, int narg, char **arg) :
-  AtomVec(lmp, narg, arg)
-{
-  molecular = 0;
-  mass_type = 1;
-  shape_type = 1;
-
-  comm_x_only = 1;
-  comm_f_only = 0;
-  size_forward = 3;
-  size_reverse = 6;
-  size_border = 6;
-  size_velocity = 6;
-  size_data_atom = 5;
-  size_data_vel = 7;
-  xcol_data = 3;
-
-  atom->omega_flag = atom->torque_flag = 1;
-}
-
-/* ----------------------------------------------------------------------
-   grow atom arrays
-   n = 0 grows arrays by DELTA
-   n > 0 allocates arrays to size n
-------------------------------------------------------------------------- */
-
-void AtomVecColloid::grow(int n)
-{
-  if (n == 0) nmax += DELTA;
-  else nmax = n;
-  atom->nmax = nmax;
-  if (nmax < 0 || nmax > MAXSMALLINT)
-    error->one("Per-processor system is too big");
-
-  tag = memory->grow(atom->tag,nmax,"atom:tag");
-  type = memory->grow(atom->type,nmax,"atom:type");
-  mask = memory->grow(atom->mask,nmax,"atom:mask");
-  image = memory->grow(atom->image,nmax,"atom:image");
-  x = memory->grow(atom->x,nmax,3,"atom:x");
-  v = memory->grow(atom->v,nmax,3,"atom:v");
-  f = memory->grow(atom->f,nmax,3,"atom:f");
-
-  omega = memory->grow(atom->omega,nmax,3,"atom:omega");
-  torque = memory->grow(atom->torque,nmax,3,"atom:torque");
-
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
-}
-
-/* ----------------------------------------------------------------------
-   reset local array ptrs
-------------------------------------------------------------------------- */
-
-void AtomVecColloid::grow_reset()
-{
-  tag = atom->tag; type = atom->type;
-  mask = atom->mask; image = atom->image;
-  x = atom->x; v = atom->v; f = atom->f;
-  omega = atom->omega; torque = atom->torque;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecColloid::copy(int i, int j)
-{
-  tag[j] = tag[i];
-  type[j] = type[i];
-  mask[j] = mask[i];
-  image[j] = image[i];
-  x[j][0] = x[i][0];
-  x[j][1] = x[i][1];
-  x[j][2] = x[i][2];
-  v[j][0] = v[i][0];
-  v[j][1] = v[i][1];
-  v[j][2] = v[i][2];
-
-  omega[j][0] = omega[i][0];
-  omega[j][1] = omega[i][1];
-  omega[j][2] = omega[i][2];
-
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_comm(int n, int *list, double *buf,
-			     int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
-      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
-      dz = pbc[2]*domain->zprd;
-    }
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0] + dx;
-      buf[m++] = x[j][1] + dy;
-      buf[m++] = x[j][2] + dz;
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_comm_vel(int n, int *list, double *buf,
-				  int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz,dvx,dvy,dvz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = v[j][0];
-      buf[m++] = v[j][1];
-      buf[m++] = v[j][2];
-      buf[m++] = omega[j][0];
-      buf[m++] = omega[j][1];
-      buf[m++] = omega[j][2];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
-      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
-      dz = pbc[2]*domain->zprd;
-    }
-    if (!deform_vremap) {
-      for (i = 0; i < n; i++) {
-	j = list[i];
-	buf[m++] = x[j][0] + dx;
-	buf[m++] = x[j][1] + dy;
-	buf[m++] = x[j][2] + dz;
-	buf[m++] = v[j][0];
-	buf[m++] = v[j][1];
-	buf[m++] = v[j][2];
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
-      }
-    } else {
-      dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
-      dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
-      dvz = pbc[2]*h_rate[2];
-      for (i = 0; i < n; i++) {
-	j = list[i];
-	buf[m++] = x[j][0] + dx;
-	buf[m++] = x[j][1] + dy;
-	buf[m++] = x[j][2] + dz;
-	if (mask[i] & deform_groupbit) {
-	  buf[m++] = v[j][0] + dvx;
-	  buf[m++] = v[j][1] + dvy;
-	  buf[m++] = v[j][2] + dvz;
-	} else {
-	  buf[m++] = v[j][0];
-	  buf[m++] = v[j][1];
-	  buf[m++] = v[j][2];
-	}
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
-      }
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecColloid::unpack_comm(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-  }
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecColloid::unpack_comm_vel(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    v[i][0] = buf[m++];
-    v[i][1] = buf[m++];
-    v[i][2] = buf[m++];
-    omega[i][0] = buf[m++];
-    omega[i][1] = buf[m++];
-    omega[i][2] = buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_reverse(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    buf[m++] = f[i][0];
-    buf[m++] = f[i][1];
-    buf[m++] = f[i][2];
-    buf[m++] = torque[i][0];
-    buf[m++] = torque[i][1];
-    buf[m++] = torque[i][2];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_reverse_one(int i, double *buf)
-{
-  buf[0] = torque[i][0];
-  buf[1] = torque[i][1];
-  buf[2] = torque[i][2];
-  return 3;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecColloid::unpack_reverse(int n, int *list, double *buf)
-{
-  int i,j,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    f[j][0] += buf[m++];
-    f[j][1] += buf[m++];
-    f[j][2] += buf[m++];
-    torque[j][0] += buf[m++];
-    torque[j][1] += buf[m++];
-    torque[j][2] += buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::unpack_reverse_one(int i, double *buf)
-{
-  torque[i][0] += buf[0];
-  torque[i][1] += buf[1];
-  torque[i][2] += buf[2];
-  return 3;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_border(int n, int *list, double *buf,
-				  int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = tag[j];
-      buf[m++] = type[j];
-      buf[m++] = mask[j];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0];
-      dy = pbc[1];
-      dz = pbc[2];
-    }
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0] + dx;
-      buf[m++] = x[j][1] + dy;
-      buf[m++] = x[j][2] + dz;
-      buf[m++] = tag[j];
-      buf[m++] = type[j];
-      buf[m++] = mask[j];
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_border_vel(int n, int *list, double *buf,
-				    int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz,dvx,dvy,dvz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = tag[j];
-      buf[m++] = type[j];
-      buf[m++] = mask[j];
-      buf[m++] = v[j][0];
-      buf[m++] = v[j][1];
-      buf[m++] = v[j][2];
-      buf[m++] = omega[j][0];
-      buf[m++] = omega[j][1];
-      buf[m++] = omega[j][2];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0];
-      dy = pbc[1];
-      dz = pbc[2];
-    }
-    if (!deform_vremap) {
-      for (i = 0; i < n; i++) {
-	j = list[i];
-	buf[m++] = x[j][0] + dx;
-	buf[m++] = x[j][1] + dy;
-	buf[m++] = x[j][2] + dz;
-	buf[m++] = tag[j];
-	buf[m++] = type[j];
-	buf[m++] = mask[j];
-	buf[m++] = v[j][0];
-	buf[m++] = v[j][1];
-	buf[m++] = v[j][2];
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
-      }
-    } else {
-      dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
-      dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
-      dvz = pbc[2]*h_rate[2];
-      for (i = 0; i < n; i++) {
-	j = list[i];
-	buf[m++] = x[j][0] + dx;
-	buf[m++] = x[j][1] + dy;
-	buf[m++] = x[j][2] + dz;
-	buf[m++] = tag[j];
-	buf[m++] = type[j];
-	buf[m++] = mask[j];
-	if (mask[i] & deform_groupbit) {
-	  buf[m++] = v[j][0] + dvx;
-	  buf[m++] = v[j][1] + dvy;
-	  buf[m++] = v[j][2] + dvz;
-	} else {
-	  buf[m++] = v[j][0];
-	  buf[m++] = v[j][1];
-	  buf[m++] = v[j][2];
-	}
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
-      }
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecColloid::unpack_border(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    if (i == nmax) grow(0);
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    tag[i] = static_cast<int> (buf[m++]);
-    type[i] = static_cast<int> (buf[m++]);
-    mask[i] = static_cast<int> (buf[m++]);
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecColloid::unpack_border_vel(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    if (i == nmax) grow(0);
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    tag[i] = static_cast<int> (buf[m++]);
-    type[i] = static_cast<int> (buf[m++]);
-    mask[i] = static_cast<int> (buf[m++]);
-    v[i][0] = buf[m++];
-    v[i][1] = buf[m++];
-    v[i][2] = buf[m++];
-    omega[i][0] = buf[m++];
-    omega[i][1] = buf[m++];
-    omega[i][2] = buf[m++];
-  }
-}
-
-/* ----------------------------------------------------------------------
-   pack data for atom I for sending to another proc
-   xyz must be 1st 3 values, so comm::exchange() can test on them
-------------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_exchange(int i, double *buf)
-{
-  int m = 1;
-  buf[m++] = x[i][0];
-  buf[m++] = x[i][1];
-  buf[m++] = x[i][2];
-  buf[m++] = v[i][0];
-  buf[m++] = v[i][1];
-  buf[m++] = v[i][2];
-  buf[m++] = tag[i];
-  buf[m++] = type[i];
-  buf[m++] = mask[i];
-  buf[m++] = image[i];
-
-  buf[m++] = omega[i][0];
-  buf[m++] = omega[i][1];
-  buf[m++] = omega[i][2];
-  
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
-
-  buf[0] = m;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecColloid::unpack_exchange(double *buf)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) grow(0);
-
-  int m = 1;
-  x[nlocal][0] = buf[m++];
-  x[nlocal][1] = buf[m++];
-  x[nlocal][2] = buf[m++];
-  v[nlocal][0] = buf[m++];
-  v[nlocal][1] = buf[m++];
-  v[nlocal][2] = buf[m++];
-  tag[nlocal] = static_cast<int> (buf[m++]);
-  type[nlocal] = static_cast<int> (buf[m++]);
-  mask[nlocal] = static_cast<int> (buf[m++]);
-  image[nlocal] = static_cast<int> (buf[m++]);
-
-  omega[nlocal][0] = buf[m++];
-  omega[nlocal][1] = buf[m++];
-  omega[nlocal][2] = buf[m++];
-  
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      m += modify->fix[atom->extra_grow[iextra]]->
-	unpack_exchange(nlocal,&buf[m]);
-
-  atom->nlocal++;
-  return m;
-}
-
-/* ----------------------------------------------------------------------
-   size of restart data for all atoms owned by this proc
-   include extra data stored by fixes
-------------------------------------------------------------------------- */
-
-int AtomVecColloid::size_restart()
-{
-  int i;
-
-  int nlocal = atom->nlocal;
-  int n = 14 * nlocal;
-
-  if (atom->nextra_restart)
-    for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
-      for (i = 0; i < nlocal; i++)
-	n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
-
-  return n;
-}
-
-/* ----------------------------------------------------------------------
-   pack atom I's data for restart file including extra quantities
-   xyz must be 1st 3 values, so that read_restart can test on them
-   molecular types may be negative, but write as positive
-------------------------------------------------------------------------- */
-
-int AtomVecColloid::pack_restart(int i, double *buf)
-{
-  int m = 1;
-  buf[m++] = x[i][0];
-  buf[m++] = x[i][1];
-  buf[m++] = x[i][2];
-  buf[m++] = tag[i];
-  buf[m++] = type[i];
-  buf[m++] = mask[i];
-  buf[m++] = image[i];
-  buf[m++] = v[i][0];
-  buf[m++] = v[i][1];
-  buf[m++] = v[i][2];
-
-  buf[m++] = omega[i][0];
-  buf[m++] = omega[i][1];
-  buf[m++] = omega[i][2];
-  
-  if (atom->nextra_restart)
-    for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
-      m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
-
-  buf[0] = m;
-  return m;
-}
-
-/* ----------------------------------------------------------------------
-   unpack data for one atom from restart file including extra quantities
-------------------------------------------------------------------------- */
-
-int AtomVecColloid::unpack_restart(double *buf)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) {
-    grow(0);
-    if (atom->nextra_store) 
-      memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
-  }
-
-  int m = 1;
-  x[nlocal][0] = buf[m++];
-  x[nlocal][1] = buf[m++];
-  x[nlocal][2] = buf[m++];
-  tag[nlocal] = static_cast<int> (buf[m++]);
-  type[nlocal] = static_cast<int> (buf[m++]);
-  mask[nlocal] = static_cast<int> (buf[m++]);
-  image[nlocal] = static_cast<int> (buf[m++]);
-  v[nlocal][0] = buf[m++];
-  v[nlocal][1] = buf[m++];
-  v[nlocal][2] = buf[m++];
-
-  omega[nlocal][0] = buf[m++];
-  omega[nlocal][1] = buf[m++];
-  omega[nlocal][2] = buf[m++];
-  
-  double **extra = atom->extra;
-  if (atom->nextra_store) {
-    int size = static_cast<int> (buf[0]) - m;
-    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
-  }
-
-  atom->nlocal++;
-  return m;
-}
-
-/* ----------------------------------------------------------------------
-   create one atom of itype at coord
-   set other values to defaults
-------------------------------------------------------------------------- */
-
-void AtomVecColloid::create_atom(int itype, double *coord)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) grow(0);
-
-  tag[nlocal] = 0;
-  type[nlocal] = itype;
-  x[nlocal][0] = coord[0];
-  x[nlocal][1] = coord[1];
-  x[nlocal][2] = coord[2];
-  mask[nlocal] = 1;
-  image[nlocal] = (512 << 20) | (512 << 10) | 512;
-  v[nlocal][0] = 0.0;
-  v[nlocal][1] = 0.0;
-  v[nlocal][2] = 0.0;
-
-  omega[nlocal][0] = 0.0;
-  omega[nlocal][1] = 0.0;
-  omega[nlocal][2] = 0.0;
-
-  atom->nlocal++;
-}
-
-/* ----------------------------------------------------------------------
-   unpack one line from Atoms section of data file
-   initialize other atom quantities
-------------------------------------------------------------------------- */
-
-void AtomVecColloid::data_atom(double *coord, int imagetmp, char **values)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) grow(0);
-
-  tag[nlocal] = atoi(values[0]);
-  if (tag[nlocal] <= 0)
-    error->one("Invalid atom ID in Atoms section of data file");
-
-  type[nlocal] = atoi(values[1]);
-  if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
-    error->one("Invalid atom type in Atoms section of data file");
-
-  x[nlocal][0] = coord[0];
-  x[nlocal][1] = coord[1];
-  x[nlocal][2] = coord[2];
-
-  image[nlocal] = imagetmp;
-
-  mask[nlocal] = 1;
-  v[nlocal][0] = 0.0;
-  v[nlocal][1] = 0.0;
-  v[nlocal][2] = 0.0;
-  omega[nlocal][0] = 0.0;
-  omega[nlocal][1] = 0.0;
-  omega[nlocal][2] = 0.0;
-
-  atom->nlocal++;
-}
-
-/* ----------------------------------------------------------------------
-   unpack one line from Velocities section of data file
-------------------------------------------------------------------------- */
-
-void AtomVecColloid::data_vel(int m, char **values)
-{
-  v[m][0] = atof(values[0]);
-  v[m][1] = atof(values[1]);
-  v[m][2] = atof(values[2]);
-  omega[m][0] = atof(values[3]);
-  omega[m][1] = atof(values[4]);
-  omega[m][2] = atof(values[5]);
-}
-
-/* ----------------------------------------------------------------------
-   unpack hybrid quantities from one line in Velocities section of data file
-------------------------------------------------------------------------- */
-
-int AtomVecColloid::data_vel_hybrid(int m, char **values)
-{
-  omega[m][0] = atof(values[0]);
-  omega[m][1] = atof(values[1]);
-  omega[m][2] = atof(values[2]);
-  return 3;
-}
-
-/* ----------------------------------------------------------------------
-   return # of bytes of allocated memory 
-------------------------------------------------------------------------- */
-
-bigint AtomVecColloid::memory_usage()
-{
-  bigint bytes = 0;
-
-  if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
-  if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
-  if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
-  if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
-  if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
-  if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
-  if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
-
-  if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
-  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
-
-  return bytes;
-}

src/COLLOID/atom_vec_colloid.h

@@ -1,66 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 ATOM_CLASS
-
-AtomStyle(colloid,AtomVecColloid)
-
-#else
-
-#ifndef LMP_ATOM_VEC_COLLOID_H
-#define LMP_ATOM_VEC_COLLOID_H
-
-#include "atom_vec.h"
-
-namespace LAMMPS_NS {
-
-class AtomVecColloid : public AtomVec {
- public:
-  AtomVecColloid(class LAMMPS *, int, char **);
-  virtual ~AtomVecColloid() {}
-  void grow(int);
-  void grow_reset();
-  void copy(int, int);
-  int pack_comm(int, int *, double *, int, int *);
-  int pack_comm_vel(int, int *, double *, int, int *);
-  void unpack_comm(int, int, double *);
-  void unpack_comm_vel(int, int, double *);
-  int pack_reverse(int, int, double *);
-  int pack_reverse_one(int, double *);
-  void unpack_reverse(int, int *, double *);
-  int unpack_reverse_one(int, double *);
-  int pack_border(int, int *, double *, int, int *);
-  int pack_border_vel(int, int *, double *, int, int *);
-  void unpack_border(int, int, double *);
-  void unpack_border_vel(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 *, int, char **);
-  void data_vel(int, char **);
-  int data_vel_hybrid(int, char **);
-  bigint memory_usage();
-
- private:
-  int *tag,*type,*mask,*image;
-  double **x,**v,**f;
-  double **omega,**torque;
-};
-
-}
-
-#endif
-#endif

src/COLLOID/fix_wall_colloid.cpp

@@ -35,31 +35,19 @@ FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) :
 
 void FixWallColloid::init()
 {
-  if (!atom->avec->shape_type)
+  if (!atom->sphere_flag) 
-    error->all("Fix wall/colloid requires atom attribute shape");
+    error->all("Fix wall/colloid requires atom style sphere");
-  if (atom->radius_flag)
-    error->all("Fix wall/colloid cannot be used with atom attribute diameter");
-
-  // insure all particle shapes are spherical
-  // can be polydisperse
-
-  for (int i = 1; i <= atom->ntypes; i++)
-    if ((atom->shape[i][0] != atom->shape[i][1]) || 
-	(atom->shape[i][0] != atom->shape[i][2]) ||
-	(atom->shape[i][1] != atom->shape[i][2]))
-      error->all("Fix wall/colloid requires spherical particles");
 
   // insure all particles in group are extended particles
 
-  double **shape = atom->shape;
+  double *radius = atom->radius;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   int flag = 0;
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit)
-      if (atom->shape[type[i]][0] == 0.0) flag = 1;
+      if (radius[i] == 0.0) flag = 1;
 
   int flagall;
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
@@ -95,9 +83,8 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
 
   double **x = atom->x;
   double **f = atom->f;
-  double **shape = atom->shape;
+  double *radius = atom->radius;
   int *mask = atom->mask;
-  int *type = atom->type;
   int nlocal = atom->nlocal;
 
   int dim = which  / 2;
@@ -111,7 +98,7 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
       if (side < 0) delta = x[i][dim] - coord;
       else delta = coord - x[i][dim];
       if (delta >= cutoff[m]) continue;
-      rad = shape[type[i]][0];
+      rad = radius[i];
       if (rad >= delta) {
 	onflag = 1;
 	continue;

src/COLLOID/pair_colloid.cpp

@@ -23,6 +23,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
+#include "neighbor.h"
 #include "neigh_list.h"
 #include "memory.h"
 #include "error.h"

src/COLLOID/pair_lubricate.cpp

@@ -85,9 +85,9 @@ void PairLubricate::compute(int eflag, int vflag)
   double **v = atom->v;
   double **f = atom->f;
   double **omega = atom->omega;
-  double **angmom = atom->angmom;
   double **torque = atom->torque;
-  double **shape = atom->shape;
+  double *radius = atom->radius;
+  double *rmass = atom->rmass;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
@@ -111,7 +111,7 @@ void PairLubricate::compute(int eflag, int vflag)
     ytmp = x[i][1];
     ztmp = x[i][2];
     itype = type[i];
-    radi = shape[itype][0];
+    radi = radius[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
@@ -151,19 +151,10 @@ void PairLubricate::compute(int eflag, int vflag)
         vt3 = vr3 - vn3;
 
         // additive rotational velocity = omega_1 + omega_2
-	// use omega directly if it exists, else angmom
-	// angular momentum = I*omega = 2/5 * M*R^2 * omega
 
-	if (omega_flag) {
+	w1 = omega[i][0] + omega[j][0];
-	  w1 = omega[i][0] + omega[j][0];
+	w2 = omega[i][1] + omega[j][1];
-	  w2 = omega[i][1] + omega[j][1];
+	w3 = omega[i][2] + omega[j][2];
-	  w3 = omega[i][2] + omega[j][2];
-	} else {
-	  inv_inertia = 1.0 / (0.4*atom->mass[itype]*radi*radi);
-	  w1 = inv_inertia * (angmom[i][0] + angmom[j][0]);
-	  w2 = inv_inertia * (angmom[i][1] + angmom[j][1]);
-	  w3 = inv_inertia * (angmom[i][2] + angmom[j][2]);
-	}
 
         // relative velocities n X P . (v1-v2) = n X (I-nn) . (v1-v2)
 
@@ -185,15 +176,9 @@ void PairLubricate::compute(int eflag, int vflag)
 
         // N.(w1-w2) and P.(w1-w2)
 
-	if (omega_flag) {
+	wr1 = omega[i][0] - omega[j][0];
-	  wr1 = omega[i][0] - omega[j][0];
+	wr2 = omega[i][1] - omega[j][1];
-	  wr2 = omega[i][1] - omega[j][1];
+	wr3 = omega[i][2] - omega[j][2];
-	  wr3 = omega[i][2] - omega[j][2];
-	} else {
-	  wr1 = inv_inertia * (angmom[i][0] - angmom[j][0]);
-	  wr2 = inv_inertia * (angmom[i][1] - angmom[j][1]);
-	  wr3 = inv_inertia * (angmom[i][2] - angmom[j][2]);
-	}
  
 	wnnr = wr1*delx + wr2*dely + wr3*delz;
 	wn1 = delx*wnnr / rsq;
@@ -393,29 +378,24 @@ void PairLubricate::coeff(int narg, char **arg)
 
 void PairLubricate::init_style()
 {
+  if (!atom->sphere_flag)
+    error->all("Pair lubricate requires atom style sphere");
   if (comm->ghost_velocity == 0)
     error->all("Pair lubricate requires ghost atoms store velocity");
-  if (!atom->torque_flag || !atom->avec->shape_type)
-    error->all("Pair lubricate requires atom attributes torque and shape");
-  if (atom->radius_flag || atom->rmass_flag)
-    error->all("Pair lubricate cannot be used with atom attributes "
-	       "diameter or rmass");
 
-  if (atom->omega_flag) omega_flag = 1;
-  else if (atom->angmom_flag) omega_flag = 0;
-  else error->all("Pair lubricate requires atom attribute omega or angmom");
-
-  // insure all particle shapes are finite-size, spherical, and monodisperse
-  // for pair hybrid, should limit test to types using the pair style
-
-  double radi = atom->shape[1][0];
-  if (radi == 0.0) error->all("Pair lubricate requires extended particles");
-  for (int i = 1; i <= atom->ntypes; i++)
-    if (atom->shape[i][0] != radi || atom->shape[i][0] != radi || 
-	atom->shape[i][0] != radi)
-      error->all("Pair lubricate requires spherical, mono-disperse particles");
-  
   int irequest = neighbor->request(this);
+
+  // require that atom radii are identical within each type require
+  // monodisperse system with same radii for all types
+
+  double rad,radtype;
+  for (int i = 1; i <= atom->ntypes; i++) {
+    if (!atom->radius_consistency(i,radtype))
+      error->all("Pair colloid requires atoms with same type have same radius");
+    if (i > 1 && radtype != rad)
+      error->all("Pair colloid requires mono-disperse particles");
+    rad = radtype;
+  }
 }
 
 /* ----------------------------------------------------------------------

src/COLLOID/pair_lubricate.h

@@ -43,7 +43,7 @@ class PairLubricate : public Pair {
   double cut_inner_global,cut_global;
   double t_target,mu;
   int flag1,flag2,flag3,flag4;
-  int seed,omega_flag;
+  int seed;
   double **cut_inner,**cut;
 
   class RanMars *random;

src/COLLOID/pair_yukawa_colloid.cpp

@@ -48,7 +48,7 @@ void PairYukawaColloid::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double **shape = atom->shape;
+  double *radius = atom->radius;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
@@ -67,7 +67,7 @@ void PairYukawaColloid::compute(int eflag, int vflag)
     ytmp = x[i][1];
     ztmp = x[i][2];
     itype = type[i];
-    radi = shape[itype][0];
+    radi = radius[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
@@ -81,7 +81,7 @@ void PairYukawaColloid::compute(int eflag, int vflag)
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-      radj = shape[jtype][0];
+      radj = radius[j];
     
       if (rsq < cutsq[itype][jtype]) {
 	r2inv = 1.0/rsq;
@@ -121,22 +121,17 @@ void PairYukawaColloid::compute(int eflag, int vflag)
 
 void PairYukawaColloid::init_style()
 {
-  if (!atom->avec->shape_type)
+  if (!atom->sphere_flag)
-    error->all("Pair yukawa/colloid requires atom attribute shape");
+    error->all("Pair yukawa/colloid requires atom style sphere");
-  if (atom->radius_flag)
-    error->all("Pair yukawa/colloid cannot be used with "
-	       "atom attribute diameter");
   
-  // insure all particle shapes are spherical
+  int irequest = neighbor->request(this);
-  // can be point particles or polydisperse
+
+  // require that atom radii are identical within each type
 
   for (int i = 1; i <= atom->ntypes; i++)
-    if ((atom->shape[i][0] != atom->shape[i][1]) || 
+    if (!atom->radius_consistency(i,rad[i]))
-	(atom->shape[i][0] != atom->shape[i][2]) ||
+      error->all("Pair yukawa/colloid requires atoms with same type "
-	(atom->shape[i][1] != atom->shape[i][2]))
+		 "have same radius");
-      error->all("Pair yukawa/colloid requires spherical particles");
-
-  int irequest = neighbor->request(this);
 }
 
 /* ----------------------------------------------------------------------
@@ -151,9 +146,7 @@ double PairYukawaColloid::init_one(int i, int j)
   }
 
   if (offset_flag) {
-    double radi = atom->shape[i][0];
+    double screening = exp(-kappa * (cut[i][j] - (rad[i]+rad[j])));
-    double radj = atom->shape[j][0];
-    double screening = exp(-kappa * (cut[i][j] - (radi+radj)));
     offset[i][j] = a[i][j]/kappa * screening;
   } else offset[i][j] = 0.0;
 
@@ -170,16 +163,12 @@ double PairYukawaColloid::single(int i, int j, int itype, int jtype,
 				 double factor_coul, double factor_lj,
 				 double &fforce)
 {
-  double r2inv,r,rinv,screening,forceyukawa,phi,radi,radj;
+  double r2inv,r,rinv,screening,forceyukawa,phi;
-
+  
-  int *type = atom->type;
-  radi = atom->shape[itype][0];
-  radj = atom->shape[jtype][0];
-
   r2inv = 1.0/rsq;
   r = sqrt(rsq);
   rinv = 1.0/r;
-  screening = exp(-kappa*(r-(radi+radj)));
+  screening = exp(-kappa*(r-(rad[itype]+rad[jtype])));
   forceyukawa = a[itype][jtype] * screening;
   fforce = factor_coul*forceyukawa * rinv;
 

src/DIPOLE/atom_vec_dipole.cpp

@@ -12,6 +12,7 @@
 ------------------------------------------------------------------------- */
 
 #include "lmptype.h"
+#include "math.h"
 #include "stdlib.h"
 #include "atom_vec_dipole.h"
 #include "atom.h"
@@ -32,19 +33,19 @@ AtomVecDipole::AtomVecDipole(LAMMPS *lmp, int narg, char **arg) :
 {
   molecular = 0;
   mass_type = 1;
-  shape_type = 1;
-  dipole_type = 1;
 
-  comm_x_only = comm_f_only = 0;
+  comm_x_only = 0;
+  comm_f_only = 1;
   size_forward = 6;
-  size_reverse = 6;
+  size_reverse = 3;
-  size_border = 10;
+  size_border = 11;
-  size_velocity = 6;
+  size_velocity = 3;
   size_data_atom = 9;
   size_data_vel = 7;
   xcol_data = 4;
 
-  atom->q_flag = atom->mu_flag = atom->omega_flag = atom->torque_flag = 1;
+  atom->dipole_flag = 1;
+  atom->q_flag = atom->mu_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
@@ -70,9 +71,7 @@ void AtomVecDipole::grow(int n)
   f = memory->grow(atom->f,nmax,3,"atom:f");
 
   q = memory->grow(atom->q,nmax,"atom:q");
-  mu = memory->grow(atom->mu,nmax,3,"atom:mu");
+  mu = memory->grow(atom->mu,nmax,4,"atom:mu");
-  omega = memory->grow(atom->omega,nmax,3,"atom:omega");
-  torque = memory->grow(atom->torque,nmax,3,"atom:torque");
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -89,7 +88,6 @@ void AtomVecDipole::grow_reset()
   mask = atom->mask; image = atom->image;
   x = atom->x; v = atom->v; f = atom->f;
   q = atom->q; mu = atom->mu;
-  omega = atom->omega; torque = atom->torque;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -111,9 +109,7 @@ void AtomVecDipole::copy(int i, int j)
   mu[j][0] = mu[i][0];
   mu[j][1] = mu[i][1];
   mu[j][2] = mu[i][2];
-  omega[j][0] = omega[i][0];
+  mu[j][3] = mu[i][3];
-  omega[j][1] = omega[i][1];
-  omega[j][2] = omega[i][2];
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -183,9 +179,6 @@ int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf,
       buf[m++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
-      buf[m++] = omega[j][0];
-      buf[m++] = omega[j][1];
-      buf[m++] = omega[j][2];
     }
   } else {
     if (domain->triclinic == 0) {
@@ -209,9 +202,6 @@ int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf,
 	buf[m++] = v[j][0];
 	buf[m++] = v[j][1];
 	buf[m++] = v[j][2];
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
       }
     } else {
       dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
@@ -234,9 +224,6 @@ int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf,
 	  buf[m++] = v[j][1];
 	  buf[m++] = v[j][2];
 	}
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
       }
     }
   }
@@ -289,9 +276,6 @@ void AtomVecDipole::unpack_comm_vel(int n, int first, double *buf)
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
-    omega[i][0] = buf[m++];
-    omega[i][1] = buf[m++];
-    omega[i][2] = buf[m++];
   }
 }
 
@@ -317,25 +301,12 @@ int AtomVecDipole::pack_reverse(int n, int first, double *buf)
     buf[m++] = f[i][0];
     buf[m++] = f[i][1];
     buf[m++] = f[i][2];
-    buf[m++] = torque[i][0];
-    buf[m++] = torque[i][1];
-    buf[m++] = torque[i][2];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecDipole::pack_reverse_one(int i, double *buf)
-{
-  buf[0] = torque[i][0];
-  buf[1] = torque[i][1];
-  buf[2] = torque[i][2];
-  return 3;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void AtomVecDipole::unpack_reverse(int n, int *list, double *buf)
 {
   int i,j,m;
@@ -346,24 +317,11 @@ void AtomVecDipole::unpack_reverse(int n, int *list, double *buf)
     f[j][0] += buf[m++];
     f[j][1] += buf[m++];
     f[j][2] += buf[m++];
-    torque[j][0] += buf[m++];
-    torque[j][1] += buf[m++];
-    torque[j][2] += buf[m++];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecDipole::unpack_reverse_one(int i, double *buf)
-{
-  torque[i][0] += buf[0];
-  torque[i][1] += buf[1];
-  torque[i][2] += buf[2];
-  return 3;
-}
-
-/* ---------------------------------------------------------------------- */
-
 int AtomVecDipole::pack_border(int n, int *list, double *buf,
 			       int pbc_flag, int *pbc)
 {
@@ -384,6 +342,7 @@ int AtomVecDipole::pack_border(int n, int *list, double *buf,
       buf[m++] = mu[j][0];
       buf[m++] = mu[j][1];
       buf[m++] = mu[j][2];
+      buf[m++] = mu[j][3];
     }
   } else {
     if (domain->triclinic == 0) {
@@ -407,6 +366,7 @@ int AtomVecDipole::pack_border(int n, int *list, double *buf,
       buf[m++] = mu[j][0];
       buf[m++] = mu[j][1];
       buf[m++] = mu[j][2];
+      buf[m++] = mu[j][3];
     }
   }
   return m;
@@ -434,12 +394,10 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf,
       buf[m++] = mu[j][0];
       buf[m++] = mu[j][1];
       buf[m++] = mu[j][2];
+      buf[m++] = mu[j][3];
       buf[m++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
-      buf[m++] = omega[j][0];
-      buf[m++] = omega[j][1];
-      buf[m++] = omega[j][2];
     }
   } else {
     if (domain->triclinic == 0) {
@@ -464,12 +422,10 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf,
 	buf[m++] = mu[j][0];
 	buf[m++] = mu[j][1];
 	buf[m++] = mu[j][2];
+	buf[m++] = mu[j][3];
 	buf[m++] = v[j][0];
 	buf[m++] = v[j][1];
 	buf[m++] = v[j][2];
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
       }
     } else {
       dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
@@ -487,6 +443,7 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf,
 	buf[m++] = mu[j][0];
 	buf[m++] = mu[j][1];
 	buf[m++] = mu[j][2];
+	buf[m++] = mu[j][3];
 	if (mask[i] & deform_groupbit) {
 	  buf[m++] = v[j][0] + dvx;
 	  buf[m++] = v[j][1] + dvy;
@@ -496,9 +453,6 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf,
 	  buf[m++] = v[j][1];
 	  buf[m++] = v[j][2];
 	}
-	buf[m++] = omega[j][0];
-	buf[m++] = omega[j][1];
-	buf[m++] = omega[j][2];
       }
     }
   }
@@ -513,7 +467,8 @@ int AtomVecDipole::pack_border_one(int i, double *buf)
   buf[1] = mu[i][0];
   buf[2] = mu[i][1];
   buf[3] = mu[i][2];
-  return 4;
+  buf[4] = mu[i][3];
+  return 5;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -536,6 +491,7 @@ void AtomVecDipole::unpack_border(int n, int first, double *buf)
     mu[i][0] = buf[m++];
     mu[i][1] = buf[m++];
     mu[i][2] = buf[m++];
+    mu[i][3] = buf[m++];
   }
 }
 
@@ -559,12 +515,10 @@ void AtomVecDipole::unpack_border_vel(int n, int first, double *buf)
     mu[i][0] = buf[m++];
     mu[i][1] = buf[m++];
     mu[i][2] = buf[m++];
+    mu[i][3] = buf[m++];
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
-    omega[i][0] = buf[m++];
-    omega[i][1] = buf[m++];
-    omega[i][2] = buf[m++];
   }
 }
 
@@ -576,7 +530,8 @@ int AtomVecDipole::unpack_border_one(int i, double *buf)
   mu[i][0] = buf[1];
   mu[i][1] = buf[2];
   mu[i][2] = buf[3];
-  return 4;
+  mu[i][3] = buf[4];
+  return 5;
 }
 
 /* ----------------------------------------------------------------------
@@ -602,9 +557,7 @@ int AtomVecDipole::pack_exchange(int i, double *buf)
   buf[m++] = mu[i][0];
   buf[m++] = mu[i][1];
   buf[m++] = mu[i][2];
-  buf[m++] = omega[i][0];
+  buf[m++] = mu[i][3];
-  buf[m++] = omega[i][1];
-  buf[m++] = omega[i][2];
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -637,9 +590,7 @@ int AtomVecDipole::unpack_exchange(double *buf)
   mu[nlocal][0] = buf[m++];
   mu[nlocal][1] = buf[m++];
   mu[nlocal][2] = buf[m++];
-  omega[nlocal][0] = buf[m++];
+  mu[nlocal][3] = buf[m++];
-  omega[nlocal][1] = buf[m++];
-  omega[nlocal][2] = buf[m++];
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -660,7 +611,7 @@ int AtomVecDipole::size_restart()
   int i;
 
   int nlocal = atom->nlocal;
-  int n = 18 * nlocal;
+  int n = 15 * nlocal;
 
   if (atom->nextra_restart)
     for (int iextra = 0; iextra < atom->nextra_restart; iextra++) 
@@ -694,9 +645,6 @@ int AtomVecDipole::pack_restart(int i, double *buf)
   buf[m++] = mu[i][0];
   buf[m++] = mu[i][1];
   buf[m++] = mu[i][2];
-  buf[m++] = omega[i][0];
-  buf[m++] = omega[i][1];
-  buf[m++] = omega[i][2];
 
   if (atom->nextra_restart)
     for (int iextra = 0; iextra < atom->nextra_restart; iextra++) 
@@ -735,9 +683,6 @@ int AtomVecDipole::unpack_restart(double *buf)
   mu[nlocal][0] = buf[m++];
   mu[nlocal][1] = buf[m++];
   mu[nlocal][2] = buf[m++];
-  omega[nlocal][0] = buf[m++];
-  omega[nlocal][1] = buf[m++];
-  omega[nlocal][2] = buf[m++];
 
   double **extra = atom->extra;
   if (atom->nextra_store) {
@@ -774,9 +719,7 @@ void AtomVecDipole::create_atom(int itype, double *coord)
   mu[nlocal][0] = 0.0;
   mu[nlocal][1] = 0.0;
   mu[nlocal][2] = 0.0;
-  omega[nlocal][0] = 0.0;
+  mu[nlocal][3] = 0.0;
-  omega[nlocal][1] = 0.0;
-  omega[nlocal][2] = 0.0;
 
   atom->nlocal++;
 }
@@ -808,6 +751,9 @@ void AtomVecDipole::data_atom(double *coord, int imagetmp, char **values)
   mu[nlocal][0] = atof(values[6]);
   mu[nlocal][1] = atof(values[7]);
   mu[nlocal][2] = atof(values[8]);
+  mu[nlocal][3] = sqrt(mu[nlocal][0]*mu[nlocal][0] + 
+		       mu[nlocal][1]*mu[nlocal][1] + 
+		       mu[nlocal][2]*mu[nlocal][2]);
 
   image[nlocal] = imagetmp;
 
@@ -815,9 +761,6 @@ void AtomVecDipole::data_atom(double *coord, int imagetmp, char **values)
   v[nlocal][0] = 0.0;
   v[nlocal][1] = 0.0;
   v[nlocal][2] = 0.0;
-  omega[nlocal][0] = 0.0;
-  omega[nlocal][1] = 0.0;
-  omega[nlocal][2] = 0.0;
 
   atom->nlocal++;
 }
@@ -833,36 +776,12 @@ int AtomVecDipole::data_atom_hybrid(int nlocal, char **values)
   mu[nlocal][0] = atof(values[1]);
   mu[nlocal][1] = atof(values[2]);
   mu[nlocal][2] = atof(values[3]);
-
+  mu[nlocal][3] = sqrt(mu[nlocal][0]*mu[nlocal][0] + 
+		       mu[nlocal][1]*mu[nlocal][1] + 
+		       mu[nlocal][2]*mu[nlocal][2]);
   return 4;
 }
 
-/* ----------------------------------------------------------------------
-   unpack one line from Velocities section of data file
-------------------------------------------------------------------------- */
-
-void AtomVecDipole::data_vel(int m, char **values)
-{
-  v[m][0] = atof(values[0]);
-  v[m][1] = atof(values[1]);
-  v[m][2] = atof(values[2]);
-  omega[m][0] = atof(values[3]);
-  omega[m][1] = atof(values[4]);
-  omega[m][2] = atof(values[5]);
-}
-
-/* ----------------------------------------------------------------------
-   unpack hybrid quantities from one line in Velocities section of data file
-------------------------------------------------------------------------- */
-
-int AtomVecDipole::data_vel_hybrid(int m, char **values)
-{
-  omega[m][0] = atof(values[0]);
-  omega[m][1] = atof(values[1]);
-  omega[m][2] = atof(values[2]);
-  return 3;
-}
-
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory 
 ------------------------------------------------------------------------- */
@@ -880,9 +799,7 @@ bigint AtomVecDipole::memory_usage()
   if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
 
   if (atom->memcheck("q")) bytes += memory->usage(q,nmax);
-  if (atom->memcheck("mu")) bytes += memory->usage(mu,nmax,3);
+  if (atom->memcheck("mu")) bytes += memory->usage(mu,nmax,4);
-  if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
-  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
 
   return bytes;
 }

src/DIPOLE/atom_vec_dipole.h

@@ -37,9 +37,7 @@ class AtomVecDipole : public AtomVec {
   void unpack_comm_vel(int, int, double *);
   int unpack_comm_one(int, double *);
   int pack_reverse(int, int, double *);
-  int pack_reverse_one(int, double *);
   void unpack_reverse(int, int *, double *);
-  int unpack_reverse_one(int, double *);
   int pack_border(int, int *, double *, int, int *);
   int pack_border_vel(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
@@ -54,8 +52,6 @@ class AtomVecDipole : public AtomVec {
   void create_atom(int, double *);
   void data_atom(double *, int, char **);
   int data_atom_hybrid(int, char **);
-  void data_vel(int, char **);
-  int data_vel_hybrid(int, char **);
   bigint memory_usage();
 
  private:

src/DIPOLE/pair_dipole_cut.cpp

@@ -79,7 +79,6 @@ void PairDipoleCut::compute(int eflag, int vflag)
   double **mu = atom->mu;
   double **torque = atom->torque;
   int *type = atom->type;
-  double *dipole = atom->dipole;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
   double *special_lj = force->special_lj;
@@ -137,7 +136,7 @@ void PairDipoleCut::compute(int eflag, int vflag)
 	    forcecoulz += pre1*delz;
 	  }
 
-	  if (dipole[itype] > 0.0 && dipole[jtype] > 0.0) { 
+	  if (mu[i][3] > 0.0 && mu[j][3] > 0.0) { 
             r3inv = r2inv*rinv;
             r5inv = r3inv*r2inv;
 	    r7inv = r5inv*r2inv;
@@ -167,7 +166,7 @@ void PairDipoleCut::compute(int eflag, int vflag)
 	    tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx);
 	  }
 
-	  if (dipole[itype] > 0.0 && q[j] != 0.0) { 
+	  if (mu[i][3] > 0.0 && q[j] != 0.0) { 
             r3inv = r2inv*rinv;
             r5inv = r3inv*r2inv;
             pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
@@ -182,7 +181,7 @@ void PairDipoleCut::compute(int eflag, int vflag)
 	    tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx);
 	  }
 
-	  if (dipole[jtype] > 0.0 && qtmp != 0.0) { 
+	  if (mu[j][3] > 0.0 && qtmp != 0.0) { 
             r3inv = r2inv*rinv;
             r5inv = r3inv*r2inv;
             pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
@@ -234,11 +233,11 @@ void PairDipoleCut::compute(int eflag, int vflag)
 	if (eflag) {
 	  if (rsq < cut_coulsq[itype][jtype]) {
 	    ecoul = qtmp*q[j]*rinv;
-	    if (dipole[itype] > 0.0 && dipole[jtype] > 0.0)
+	    if (mu[i][3] > 0.0 && mu[j][3] > 0.0)
 	      ecoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr;
-	    if (dipole[itype] > 0.0 && q[j] != 0.0) 
+	    if (mu[i][3] > 0.0 && q[j] != 0.0) 
 	      ecoul += -q[j]*r3inv*pidotr;
-	    if (dipole[jtype] > 0.0 && qtmp != 0.0)
+	    if (mu[j][3] > 0.0 && qtmp != 0.0)
 	      ecoul += qtmp*r3inv*pjdotr;
 	    ecoul *= factor_coul*qqrd2e;
 	  } else ecoul = 0.0;
@@ -357,10 +356,8 @@ void PairDipoleCut::coeff(int narg, char **arg)
 
 void PairDipoleCut::init_style()
 {
-  if (!atom->q_flag || !atom->mu_flag || 
+  if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag)
-      !atom->torque_flag || atom->dipole == NULL)
+    error->all("Pair dipole/cut requires atom attributes q, mu, torque");
-    error->all("Pair dipole/cut requires atom attributes "
-	       "q, mu, torque, dipole");
 
   int irequest = neighbor->request(this);
 }

src/GPU/pair_gayberne_gpu.cpp

@@ -134,19 +134,21 @@ void PairGayBerneGPU::init_style()
 {
   if (force->pair_match("gpu",0) == NULL)
     error->all("Cannot use pair hybrid with multiple GPU pair styles");
-  if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type)
+  if (!atom->ellipsoid_flag)
-    error->all("Pair gayberne requires atom attributes quat, torque, shape");
+    error->all("Pair gayberne requires atom style ellipsoid");
-  if (atom->radius_flag)
-    error->all("Pair gayberne cannot be used with atom attribute diameter");
 
   // per-type shape precalculations
+  // require that atom shapes are identical within each type
+
   for (int i = 1; i <= atom->ntypes; i++) {
+    if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2]))
+      error->all("Pair gayberne requires atoms with same type have same shape");
     if (setwell[i]) {
-      double *one = atom->shape[i];
+      shape2[i][0] = shape1[i][0]*shape1[i][0];
-      shape[i][0] = one[0]*one[0];
+      shape2[i][1] = shape1[i][1]*shape1[i][1];
-      shape[i][1] = one[1]*one[1];
+      shape2[i][2] = shape1[i][2]*shape1[i][2];
-      shape[i][2] = one[2]*one[2];
+      lshape[i] = (shape1[i][0]*shape1[i][1]+shape1[i][2]*shape1[i][2]) * 
-      lshape[i] = (one[0]*one[1]+one[2]*one[2])*sqrt(one[0]*one[1]);
+	sqrt(shape1[i][0]*shape1[i][1]);
     }
   }
 
@@ -169,7 +171,7 @@ void PairGayBerneGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   bool init_ok = gb_gpu_init(atom->ntypes+1, gamma, upsilon, mu, 
-                             shape, well, cutsq, sigma, epsilon, lshape, form,
+                             shape1, well, cutsq, sigma, epsilon, lshape, form,
                              lj1, lj2, lj3, lj4, offset, force->special_lj, 
                              atom->nlocal, atom->nlocal+atom->nghost, 300, 
                              cell_size, gpu_mode, screen);
@@ -227,7 +229,7 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag)
       MathExtra::quat_to_mat_trans(quat[i],a1);
       MathExtra::diag_times3(well[itype],a1,temp);
       MathExtra::transpose_times3(a1,temp,b1);
-      MathExtra::diag_times3(shape[itype],a1,temp);
+      MathExtra::diag_times3(shape2[itype],a1,temp);
       MathExtra::transpose_times3(a1,temp,g1);
     }
 
@@ -271,7 +273,7 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag)
 	  MathExtra::quat_to_mat_trans(quat[j],a2);
 	  MathExtra::diag_times3(well[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,b2);
-	  MathExtra::diag_times3(shape[jtype],a2,temp);
+	  MathExtra::diag_times3(shape2[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,g2);
 	  one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor);
 	  ttor[0] = ttor[1] = ttor[2] = 0.0;
@@ -286,7 +288,7 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag)
 	  MathExtra::quat_to_mat_trans(quat[j],a2);
 	  MathExtra::diag_times3(well[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,b2);
-	  MathExtra::diag_times3(shape[jtype],a2,temp);
+	  MathExtra::diag_times3(shape2[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,g2);
 	  one_eng = gayberne_analytic(i,j,a1,a2,b1,b2,g1,g2,r12,rsq,
 				      fforce,ttor,rtor);
@@ -343,7 +345,7 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag)
       MathExtra::quat_to_mat_trans(quat[i],a1);
       MathExtra::diag_times3(well[itype],a1,temp);
       MathExtra::transpose_times3(a1,temp,b1);
-      MathExtra::diag_times3(shape[itype],a1,temp);
+      MathExtra::diag_times3(shape2[itype],a1,temp);
       MathExtra::transpose_times3(a1,temp,g1);
     }
 
@@ -389,7 +391,7 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag)
 	  MathExtra::quat_to_mat_trans(quat[j],a2);
 	  MathExtra::diag_times3(well[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,b2);
-	  MathExtra::diag_times3(shape[jtype],a2,temp);
+	  MathExtra::diag_times3(shape2[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,g2);
 	  one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor);
 	  ttor[0] = ttor[1] = ttor[2] = 0.0;
@@ -404,7 +406,7 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag)
 	  MathExtra::quat_to_mat_trans(quat[j],a2);
 	  MathExtra::diag_times3(well[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,b2);
-	  MathExtra::diag_times3(shape[jtype],a2,temp);
+	  MathExtra::diag_times3(shape2[jtype],a2,temp);
 	  MathExtra::transpose_times3(a2,temp,g2);
 	  one_eng = gayberne_analytic(i,j,a1,a2,b1,b2,g1,g2,r12,rsq,
 				      fforce,ttor,rtor);

src/GRANULAR/Install.sh

@@ -2,7 +2,6 @@
 
 if (test $1 = 1) then
 
-  cp atom_vec_granular.cpp ..
   cp fix_freeze.cpp ..
   cp fix_pour.cpp ..
   cp fix_wall_gran.cpp ..
@@ -10,7 +9,6 @@ if (test $1 = 1) then
   cp pair_gran_hooke.cpp ..
   cp pair_gran_hooke_history.cpp ..
 
-  cp atom_vec_granular.h ..
   cp fix_freeze.h ..
   cp fix_pour.h ..
   cp fix_wall_gran.h ..
@@ -20,7 +18,6 @@ if (test $1 = 1) then
 
 elif (test $1 = 0) then
 
-  rm ../atom_vec_granular.cpp
   rm ../fix_freeze.cpp
   rm ../fix_pour.cpp
   rm ../fix_wall_gran.cpp
@@ -28,7 +25,6 @@ elif (test $1 = 0) then
   rm ../pair_gran_hooke.cpp
   rm ../pair_gran_hooke_history.cpp
 
-  rm ../atom_vec_granular.h
   rm ../fix_freeze.h
   rm ../fix_pour.h
   rm ../fix_wall_gran.h

src/GRANULAR/fix_pour.cpp

@@ -472,7 +472,6 @@ void FixPour::pre_exchange()
       m = atom->nlocal - 1;
       atom->type[m] = ntype;
       atom->radius[m] = radtmp;
-      atom->density[m] = denstmp;
       atom->rmass[m] = 4.0*PI/3.0 * radtmp*radtmp*radtmp * denstmp;
       atom->mask[m] = 1 | groupbit;
       atom->v[m][0] = vxtmp;

src/GRANULAR/fix_wall_gran.cpp

@@ -46,8 +46,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 10) error->all("Illegal fix wall/gran command");
 
-  if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)
+  if (!atom->sphere_flag)
-    error->all("Fix wall/gran requires atom attributes radius, omega, torque");
+    error->all("Fix wall/gran requires atom style sphere");
 
   restart_peratom = 1;
   create_attribute = 1;

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -364,8 +364,8 @@ void PairGranHookeHistory::init_style()
 
   // error and warning checks
 
-  if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)
+  if (!atom->sphere_flag)
-    error->all("Pair granular requires atom attributes radius, omega, torque");
+    error->all("Pair granular requires atom style sphere");
   if (comm->ghost_velocity == 0)
     error->all("Pair granular requires ghost atoms store velocity");
 

src/PERI/atom_vec_peri.cpp

@@ -47,7 +47,8 @@ AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) :
   size_data_vel = 4;
   xcol_data = 5;
 
-  atom->vfrac_flag = atom->density_flag = atom->rmass_flag = 1;
+  atom->peri_flag = 1;
+  atom->vfrac_flag = atom->rmass_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
@@ -73,7 +74,6 @@ void AtomVecPeri::grow(int n)
   f = memory->grow(atom->f,nmax,3,"atom:f");
 
   vfrac = memory->grow(atom->vfrac,nmax,"atom:vfrac");
-  density = memory->grow(atom->density,nmax,"atom:density");
   rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
   s0 = memory->grow(atom->s0,nmax,"atom:s0");
   x0 = memory->grow(atom->x0,nmax,3,"atom:x0");
@@ -92,7 +92,7 @@ void AtomVecPeri::grow_reset()
   tag = atom->tag; type = atom->type;
   mask = atom->mask; image = atom->image;
   x = atom->x; v = atom->v; f = atom->f;
-  vfrac = atom->vfrac; density = atom->density; rmass = atom->rmass;
+  vfrac = atom->vfrac; rmass = atom->rmass;
   s0 = atom->s0; x0 = atom->x0;
 }
 
@@ -112,7 +112,6 @@ void AtomVecPeri::copy(int i, int j)
   v[j][2] = v[i][2];
 
   vfrac[j] = vfrac[i];
-  density[j] = density[i];
   rmass[j] = rmass[i];
   s0[j] = s0[i];
   x0[j][0] = x0[i][0];
@@ -546,7 +545,6 @@ int AtomVecPeri::pack_exchange(int i, double *buf)
   buf[m++] = image[i];
 
   buf[m++] = vfrac[i];
-  buf[m++] = density[i];
   buf[m++] = rmass[i];
   buf[m++] = s0[i];
   buf[m++] = x0[i][0];
@@ -581,7 +579,6 @@ int AtomVecPeri::unpack_exchange(double *buf)
   image[nlocal] = static_cast<int> (buf[m++]);
 
   vfrac[nlocal] = buf[m++];
-  density[nlocal] = buf[m++];
   rmass[nlocal] = buf[m++];
   s0[nlocal] = buf[m++];
   x0[nlocal][0] = buf[m++];
@@ -608,7 +605,7 @@ int AtomVecPeri::size_restart()
   int i;
  
   int nlocal = atom->nlocal;
-  int n = 18 * nlocal;
+  int n = 17 * nlocal;
  
   if (atom->nextra_restart)
     for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@@ -639,7 +636,6 @@ int AtomVecPeri::pack_restart(int i, double *buf)
   buf[m++] = v[i][2];
  
   buf[m++] = vfrac[i];
-  buf[m++] = density[i];
   buf[m++] = rmass[i];
   buf[m++] = s0[i];
   buf[m++] = x0[i][0];
@@ -680,7 +676,6 @@ int AtomVecPeri::unpack_restart(double *buf)
   v[nlocal][2] = buf[m++];
  
   vfrac[nlocal] = buf[m++];
-  density[nlocal] = buf[m++];
   rmass[nlocal] = buf[m++];
   s0[nlocal] = buf[m++];
   x0[nlocal][0] = buf[m++];
@@ -719,8 +714,7 @@ void AtomVecPeri::create_atom(int itype, double *coord)
   v[nlocal][2] = 0.0;
  
   vfrac[nlocal] = 1.0;
-  density[nlocal] = 1.0;
+  rmass[nlocal] = 1.0;
-  rmass[nlocal] = density[nlocal];
   s0[nlocal] = DBL_MAX;
   x0[nlocal][0] = coord[0];
   x0[nlocal][1] = coord[1];
@@ -748,8 +742,7 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values)
     error->one("Invalid atom type in Atoms section of data file");
 
   vfrac[nlocal] = atof(values[2]);
-  density[nlocal] = atof(values[3]);
+  rmass[nlocal] = atof(values[3]);
-  rmass[nlocal] = density[nlocal];
   if (rmass[nlocal] <= 0.0) error->one("Invalid mass value");
 
   x[nlocal][0] = coord[0];
@@ -780,8 +773,7 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values)
 int AtomVecPeri::data_atom_hybrid(int nlocal, char **values)
 {
   vfrac[nlocal] = atof(values[0]);
-  density[nlocal] = atof(values[1]);
+  rmass[nlocal] = atof(values[1]);
-  rmass[nlocal] = density[nlocal];
   if (rmass[nlocal] <= 0.0) error->one("Invalid mass value");
 
   s0[nlocal] = DBL_MAX;
@@ -809,7 +801,6 @@ bigint AtomVecPeri::memory_usage()
   if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
 
   if (atom->memcheck("vfrac")) bytes += memory->usage(vfrac,nmax);
-  if (atom->memcheck("density")) bytes += memory->usage(density,nmax);
   if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
   if (atom->memcheck("s0")) bytes += memory->usage(s0,nmax);
   if (atom->memcheck("x0")) bytes += memory->usage(x0,nmax,3);

src/SRD/fix_srd.cpp

@@ -40,8 +40,6 @@ using namespace LAMMPS_NS;
 
 enum{SLIP,NOSLIP};
 enum{SPHERE,ELLIPSOID,WALL};
-enum{SPHERE_SHAPE,SPHERE_RADIUS};
-enum{ANGULAR_OMEGA,ANGULAR_ANGMOM};
 enum{INSIDE_ERROR,INSIDE_WARN,INSIDE_IGNORE};
 enum{BIG_MOVE,SRD_MOVE,SRD_ROTATE};
 enum{CUBIC_ERROR,CUBIC_WARN};
@@ -291,10 +289,8 @@ void FixSRD::init()
   if (bigexist && comm->ghost_velocity == 0)
     error->all("Fix srd requires ghost atoms store velocity");
 
-  if (bigexist && !atom->radius_flag && !atom->avec->shape_type)
+  if (bigexist && !atom->sphere_flag && !atom->ellipsoid_flag)
-    error->all("Fix SRD requires atom attribute radius or shape");
+    error->all("Fix SRD requires atom style sphere or ellipsoid");
-  if (bigexist && !atom->angmom_flag && !atom->omega_flag)
-    error->all("Fix SRD requires atom attribute angmom or omega");
   if (bigexist && atom->angmom_flag && atom->omega_flag)
     error->all("Fix SRD cannot have both atom attributes angmom and omega");
   if (bigexist && collidestyle == NOSLIP && !atom->torque_flag)
@@ -897,15 +893,6 @@ void FixSRD::reset_velocities()
     vold[1] /= n;
     vold[2] /= n;
 
-    // if (triclinic && streamflag) {
-    // xlamda = vbin[i].xctr;
-    // vstream[0] = h_rate[0]*xlamda[0] + h_rate[5]*xlamda[1] + 
-    //	h_rate[4]*xlamda[2] + h_ratelo[0];
-    // vstream[1] = h_rate[1]*xlamda[1] + h_rate[3]*xlamda[2] + h_ratelo[1];
-    //  vstream[2] = h_rate[2]*xlamda[2] + h_ratelo[2];
-    // vnew = vstream;
-    //} else vnew = vold;
-
     vnew = vold;
 
     irandom = static_cast<int> (6.0*vbin[i].random);
@@ -2109,7 +2096,6 @@ void FixSRD::parameterize()
 
   double *radius = atom->radius;
   double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
@@ -2119,24 +2105,20 @@ void FixSRD::parameterize()
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & biggroupbit) {
-      if (radius) {
+      if (radius && radius[i] > 0.0) {
-	if (radius[i] == 0.0)
-	  error->one("Big particle in fix srd cannot be point particle");
 	maxbigdiam = MAX(maxbigdiam,2.0*radius[i]);
 	minbigdiam = MIN(minbigdiam,2.0*radius[i]);
-      } else {
+      } else if (shape && shape[i][0] > 0.0) {
-	if (shape[type[i]][0] == 0.0)
+	maxbigdiam = MAX(maxbigdiam,2.0*shape[i][0]);
-	  error->one("Big particle in fix srd cannot be point particle");
+	maxbigdiam = MAX(maxbigdiam,2.0*shape[i][1]);
-	maxbigdiam = MAX(maxbigdiam,2.0*shape[type[i]][0]);
+	maxbigdiam = MAX(maxbigdiam,2.0*shape[i][2]);
-	maxbigdiam = MAX(maxbigdiam,2.0*shape[type[i]][1]);
+	minbigdiam = MIN(minbigdiam,2.0*shape[i][0]);
-	maxbigdiam = MAX(maxbigdiam,2.0*shape[type[i]][2]);
+	minbigdiam = MIN(minbigdiam,2.0*shape[i][1]);
-	minbigdiam = MIN(minbigdiam,2.0*shape[type[i]][0]);
+	minbigdiam = MIN(minbigdiam,2.0*shape[i][2]);
-	minbigdiam = MIN(minbigdiam,2.0*shape[type[i]][1]);
+	if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
-	minbigdiam = MIN(minbigdiam,2.0*shape[type[i]][2]);
-	if (shape[type[i]][0] != shape[type[i]][1] ||
-	    shape[type[i]][0] != shape[type[i]][2])
 	  any_ellipsoids = 1;
-      }
+      } else 
+	error->one("Big particle in fix srd cannot be point particle");
     }
 
   double tmp = maxbigdiam;
@@ -2159,6 +2141,7 @@ void FixSRD::parameterize()
 
   double *rmass = atom->rmass;
   double *mass = atom->mass;
+  int *type = atom->type;
 
   int flag = 0;
   mass_srd = 0.0;
@@ -2166,8 +2149,10 @@ void FixSRD::parameterize()
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       if (rmass) {
-	if (mass_srd == 0.0) mass_srd = rmass[i];
+	if (mass_srd == 0.0) {
-	else if (rmass[i] != mass_srd) flag = 1;
+	  mass_srd = rmass[i];
+	  printf("MASS SRD %g\n",rmass[i]);
+	} else if (rmass[i] != mass_srd) flag = 1;
       } else {
 	if (mass_srd == 0.0) mass_srd = mass[type[i]];
 	else if (mass[type[i]] != mass_srd) flag = 1;
@@ -2188,6 +2173,8 @@ void FixSRD::parameterize()
     temperature_srd = force->mvv2e * 
       (lamda/dt_srd)*(lamda/dt_srd) * mass_srd/force->boltz;
 
+  printf("AAA %g %g\n",mass_srd,lamda);
+
   // vmax = maximum velocity of an SRD particle
   // dmax = maximum distance an SRD can move = 4*lamda = vmax * dt_srd
 
@@ -2204,19 +2191,18 @@ void FixSRD::parameterize()
   if (dimension == 3) {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & biggroupbit) {
-	if (radius)
+	if (radius && radius[i] > 0.0)
 	  volbig += 4.0/3.0*PI*radius[i]*radius[i]*radius[i];
-	else
+	else if (shape && shape[i][0] > 0.0)
-	  volbig += 4.0/3.0*PI *
+	  volbig += 4.0/3.0*PI * shape[i][0]*shape[i][1]*shape[i][2];
-	    shape[type[i]][0]*shape[type[i]][1]*shape[type[i]][2];
       }
   } else {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & biggroupbit) {
-	if (radius)
+	if (radius && radius[i] > 0.0)
 	  volbig += PI*radius[i]*radius[i];
-	else
+	else if (shape && shape[i][0] > 0.0)
-	  volbig += PI*shape[type[i]][0]*shape[type[i]][1];
+	  volbig += PI*shape[i][0]*shape[i][1];
       }
   }
 
@@ -2236,10 +2222,7 @@ void FixSRD::parameterize()
 
   mass_big = 0.0;
   for (int i = 0; i < nlocal; i++)
-    if (mask[i] & biggroupbit) {
+    if (mask[i] & biggroupbit) mass_big += rmass[i];
-      if (rmass) mass_big += rmass[i];
-      else mass_big += mass[type[i]];
-    }
 
   tmp = mass_big;
   MPI_Allreduce(&tmp,&mass_big,1,MPI_DOUBLE,MPI_SUM,world);
@@ -2425,40 +2408,21 @@ void FixSRD::big_static()
   double **shape = atom->shape;
   int *type = atom->type;
 
-  int omega_flag = atom->omega_flag;
-
   double skinhalf = 0.5 * neighbor->skin;
 
   for (int k = 0; k < nbig; k++) {
     i = biglist[k].index;
-    if (radius) {
+    if (radius && radius[i] > 0.0) {
       biglist[k].type = SPHERE;
-      biglist[k].typesphere = SPHERE_RADIUS;
-      biglist[k].typeangular = ANGULAR_OMEGA;
       rad = radfactor*radius[i];
       biglist[k].radius = rad;
       biglist[k].radsq = rad*rad;
       biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
-    } else if (shape[type[i]][0] == shape[type[i]][1] && 
+    } else if (shape && shape[i][0] > 0.0) {
-	       shape[type[i]][0] == shape[type[i]][2]) {
-      biglist[k].type = SPHERE;
-      biglist[k].typesphere = SPHERE_SHAPE;
-
-      // either atom->omega is defined or atom->angmom, cannot both be defined
-
-      if (omega_flag) biglist[k].typeangular = ANGULAR_OMEGA;
-      else biglist[k].typeangular = ANGULAR_ANGMOM;
-
-      rad = radfactor*shape[type[i]][0];
-      biglist[k].radius = rad;
-      biglist[k].radsq = rad*rad;
-      biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
-    } else {
       biglist[k].type = ELLIPSOID;
-      biglist[k].typeangular = ANGULAR_ANGMOM;
+      arad = radfactor*shape[i][0];
-      arad = radfactor*shape[type[i]][0];
+      brad = radfactor*shape[i][1];
-      brad = radfactor*shape[type[i]][1];
+      crad = radfactor*shape[i][2];
-      crad = radfactor*shape[type[i]][2];
       biglist[k].aradsqinv = 1.0/(arad*arad);
       biglist[k].bradsqinv = 1.0/(brad*brad);
       biglist[k].cradsqinv = 1.0/(crad*crad);
@@ -2478,44 +2442,32 @@ void FixSRD::big_static()
 
 void FixSRD::big_dynamic()
 {
-  int i,itype;
+  int i;
-  double inertia;
 
   double **omega = atom->omega;
   double **angmom = atom->angmom;
   double **quat = atom->quat;
   double **shape = atom->shape;
-  double *mass = atom->mass;
+  double *rmass = atom->rmass;
-  int *type = atom->type;
 
   for (int k = 0; k < nbig; k++) {
     i = biglist[k].index;
 
-    // ellipsoid with angmom
+    // sphere with omega
-    // calculate ex,ey,ez and omega from quaternion and angmom
-
-    if (biglist[k].type == ELLIPSOID) {
-      exyz_from_q(quat[i],biglist[k].ex,biglist[k].ey,biglist[k].ez);
-      omega_from_mq(angmom[i],biglist[k].ex,biglist[k].ey,biglist[k].ez,
-		    mass[type[i]],shape[type[i]],biglist[k].omega);
-
-    // sphere with omega and shape or radius
     // set omega from atom->omega directly
 
-    } else if (biglist[k].typeangular == ANGULAR_OMEGA) {
+    if (biglist[k].type == SPHERE) {
       biglist[k].omega[0] = omega[i][0];
       biglist[k].omega[1] = omega[i][1];
       biglist[k].omega[2] = omega[i][2];
 
-    // sphere with angmom and shape
+    // ellipsoid with angmom
-    // calculate omega from angmom
+    // calculate ex,ey,ez and omega from quaternion and angmom
 
-    } else if (biglist[k].typeangular == ANGULAR_ANGMOM) {
+    } else if (biglist[k].type == ELLIPSOID) {
-      itype = type[i];
+      exyz_from_q(quat[i],biglist[k].ex,biglist[k].ey,biglist[k].ez);
-      inertia = INERTIA * mass[itype]*shape[itype][0]*shape[itype][0];
+      omega_from_mq(angmom[i],biglist[k].ex,biglist[k].ey,biglist[k].ez,
-      biglist[k].omega[0] = angmom[i][0] / inertia;
+		    rmass[i],shape[i],biglist[k].omega);
-      biglist[k].omega[1] = angmom[i][1] / inertia;
-      biglist[k].omega[2] = angmom[i][2] / inertia;
     }
   }
 }

src/SRD/fix_srd.h

@@ -91,8 +91,6 @@ class FixSRD : public Fix {
   struct Big {
     int index;                 // local index of particle/wall
     int type;                  // SPHERE or ELLIPSOID or WALL
-    int typesphere;            // SPHERE_SHAPE or SPHERE_RADIUS
-    int typeangular;           // ANGULAR_OMEGA or ANGULAR_ANGMOM
     double radius,radsq;       // radius of sphere
     double aradsqinv;          // 3 ellipsoid radii
     double bradsqinv;

src/USER-EFF/atom_vec_electron.cpp

@@ -48,6 +48,7 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp, int narg, char **arg) :
   size_data_vel = 5;
   xcol_data = 6;
   
+  atom->electron_flag = 1;
   atom->q_flag = atom->spin_flag = atom->eradius_flag = 
     atom->ervel_flag = atom->erforce_flag = 1;
 }

src/USER-EFF/compute_ke_atom_eff.cpp

@@ -43,8 +43,8 @@ ComputeKEAtomEff::ComputeKEAtomEff(LAMMPS *lmp, int narg, char **arg) :
 
   // error check
 
-  if (!atom->ervel_flag) 
+  if (!atom->electron_flag) 
-    error->all("Compute ke/atom/eff requires atom attribute ervel");
+    error->all("Compute ke/atom/eff requires atom style electron");
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-EFF/compute_ke_eff.cpp

@@ -39,8 +39,8 @@ ComputeKEEff::ComputeKEEff(LAMMPS *lmp, int narg, char **arg) :
 
   // error check
 
-  if (!atom->ervel_flag) 
+  if (!atom->electron_flag) 
-    error->all("Compute ke/eff requires atom attribute ervel");
+    error->all("Compute ke/eff requires atom style electron");
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-EFF/compute_temp_deform_eff.cpp

@@ -43,8 +43,8 @@ ComputeTempDeformEff::ComputeTempDeformEff(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 3) error->all("Illegal compute temp/deform/eff command");
 
-  if (!atom->spin_flag || !atom->ervel_flag) 
+  if (!atom->electron_flag) 
-    error->all("Compute temp/deform/eff requires atom attributes spin, ervel");
+    error->all("Compute temp/deform/eff requires atom style electron");
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;

src/USER-EFF/compute_temp_eff.cpp

@@ -35,10 +35,8 @@ using namespace LAMMPS_NS;
 ComputeTempEff::ComputeTempEff(LAMMPS *lmp, int narg, char **arg) : 
   Compute(lmp, narg, arg)
 {
-  // error check
+  if (!atom->electron_flag) 
-
+    error->all("Compute temp/eff requires atom style electron");
-  if (!atom->spin_flag || !atom->ervel_flag) 
-    error->all("Compute temp/eff requires atom attributes spin, ervel");
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;

src/USER-EFF/compute_temp_region_eff.cpp

@@ -35,8 +35,8 @@ using namespace LAMMPS_NS;
 ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) : 
   Compute(lmp, narg, arg)
 {
-  if (!atom->spin_flag || !atom->ervel_flag) 
+  if (!atom->electron_flag) 
-    error->all("Compute temp/region/eff requires atom attributes spin, ervel");
+    error->all("Compute temp/region/eff requires atom style electron");
 
   if (narg != 4) error->all("Illegal compute temp/region/eff command");
 

src/USER-EFF/fix_nh_eff.cpp

@@ -30,10 +30,8 @@ enum{NOBIAS,BIAS};
 
 FixNHEff::FixNHEff(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg)
 {
-  if (!atom->spin_flag || !atom->eradius_flag || 
+  if (!atom->electron_flag) 
-      !atom->ervel_flag || !atom->erforce_flag) 
+    error->all("Fix nvt/nph/npt/eff requires atom style electron");
-    error->all("Fix nvt/nph/npt/eff requires atom attributes "
-	       "spin, eradius, ervel, erforce");
 }
 
 /* ----------------------------------------------------------------------

src/USER-EFF/fix_nve_eff.cpp

@@ -33,10 +33,8 @@ using namespace LAMMPS_NS;
 FixNVEEff::FixNVEEff(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (!atom->spin_flag || !atom->eradius_flag || 
+  if (!atom->electron_flag) 
-      !atom->ervel_flag || !atom->erforce_flag) 
+    error->all("Fix nve/eff requires atom style electron");
-    error->all("Fix nve/eff requires atom attributes "
-               "spin, eradius, ervel, erforce");
 
   time_integrate = 1;
 }

src/atom.cpp

@@ -70,8 +70,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   molecule = NULL;
   q = NULL;
   mu = NULL;
-  quat = omega = angmom = torque = NULL;
+  quat = omega = angmom = torque = shape = NULL;
-  radius = density = rmass = NULL;
+  radius = rmass = NULL;
   vfrac = s0 = NULL;
   x0 = NULL;
 
@@ -96,23 +96,20 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   improper_type = improper_atom1 = improper_atom2 = NULL;
   improper_atom3 = improper_atom4 = NULL;
 
-  // initialize atom array existence flags
+  // initialize atom style and array existence flags
   // customize by adding new flag
 
-  molecule_flag = 0;
+  sphere_flag = ellipsoid_flag = peri_flag = dipole_flag = electron_flag = 0;
-  q_flag = mu_flag = 0;
+
-  quat_flag = omega_flag = angmom_flag = torque_flag = 0;
+  molecule_flag = q_flag = mu_flag = 0;
-  radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
+  rmass_flag = radius_flag = omega_flag = torque_flag = 0;
-  spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
+  quat_flag = shape_flag = angmom_flag = 0;
+  vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
 
   // ntype-length arrays
 
   mass = NULL;
   mass_setflag = NULL;
-  shape = NULL;
-  shape_setflag = NULL;
-  dipole = NULL;
-  dipole_setflag = NULL;
 
   // callback lists & extra restart info
 
@@ -171,12 +168,12 @@ Atom::~Atom()
   memory->destroy(q);
   memory->destroy(mu);
   memory->destroy(quat);
+  memory->destroy(shape);
   memory->destroy(omega);
   memory->destroy(angmom);
   memory->destroy(torque);
 
   memory->destroy(radius);
-  memory->destroy(density);
   memory->destroy(rmass);
   memory->destroy(vfrac);
   memory->destroy(s0);
@@ -220,10 +217,6 @@ Atom::~Atom()
 
   delete [] mass;
   delete [] mass_setflag;
-  memory->destroy(shape);
-  delete [] shape_setflag;
-  delete [] dipole;
-  delete [] dipole_setflag;
 
   // delete extra arrays
 
@@ -256,14 +249,16 @@ void Atom::create_avec(const char *style, int narg, char **arg)
   delete [] atom_style;
   if (avec) delete avec;
 
-  // unset atom array existence flags that may have been set by old avec
+  // unset atom style and array existence flags
+  // may have been set by old avec
   // customize by adding new flag
 
-  molecule_flag = 0;
+  sphere_flag = ellipsoid_flag = peri_flag = dipole_flag = electron_flag = 0;
-  q_flag = mu_flag = 0;
+
-  quat_flag = omega_flag = angmom_flag = torque_flag = 0;
+  molecule_flag = q_flag = mu_flag = 0;
-  radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
+  rmass_flag = radius_flag = omega_flag = torque_flag = 0;
-  spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
+  quat_flag = shape_flag = angmom_flag = 0;
+  vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
 
   avec = new_avec(style,narg,arg);
   int n = strlen(style) + 1;
@@ -309,8 +304,6 @@ void Atom::init()
   // check arrays that are atom type in length
 
   check_mass();
-  check_shape();
-  check_dipole();
 
   // setup of firstgroup
 
@@ -1066,16 +1059,6 @@ void Atom::allocate_type_arrays()
     mass_setflag = new int[ntypes+1];
     for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
   }
-  if (avec->shape_type) {
-    memory->create(shape,ntypes+1,3,"atom:shape");
-    shape_setflag = new int[ntypes+1];
-    for (int itype = 1; itype <= ntypes; itype++) shape_setflag[itype] = 0;
-  }    
-  if (avec->dipole_type) {
-    dipole = new double[ntypes+1];
-    dipole_setflag = new int[ntypes+1];
-    for (int itype = 1; itype <= ntypes; itype++) dipole_setflag[itype] = 0;
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1161,161 +1144,61 @@ void Atom::check_mass()
 }
 
 /* ----------------------------------------------------------------------
-   set particle shape and flag it as set
+   check that radii of all particles of itype are the same
-   called from reading of data file
+   return 1 if true, else return 0
+   also return the radius value for that type
 ------------------------------------------------------------------------- */
 
-void Atom::set_shape(const char *str)
+int Atom::radius_consistency(int itype, double &rad)
 {
-  if (shape == NULL) error->all("Cannot set shape for this atom style");
+  double value = -1.0;
-
+  int flag = 0;
-  int itype;
+  for (int i = 0; i < nlocal; i++) {
-  double a,b,c;
+    if (type[i] != itype) continue;
-  int n = sscanf(str,"%d %lg %lg %lg",&itype,&a,&b,&c);
+    if (value < 0.0) value = radius[i];
-  if (n != 4) error->all("Invalid shape line in data file");
+    else if (value != radius[i]) flag = 1;
-
-  if (itype < 1 || itype > ntypes) error->all("Invalid type for shape set");
-
-  // store shape as radius, though specified as diameter
-
-  shape[itype][0] = 0.5*a;
-  shape[itype][1] = 0.5*b;
-  shape[itype][2] = 0.5*c;
-  shape_setflag[itype] = 1;
-
-  if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || 
-      shape[itype][2] < 0.0)
-    error->all("Invalid shape value");
-  if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || 
-      shape[itype][2] > 0.0) {
-    if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || 
-	shape[itype][2] == 0.0)
-      error->all("Invalid shape value");
   }
+
+  int flagall;
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+  if (flagall) return 0;
+
+  MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world);
+  return 1;
 }
 
 /* ----------------------------------------------------------------------
-   set one or more particle shapes and flag them as set
+   check that shape of all particles of itype are the same
-   called from reading of input script
+   return 1 if true, else return 0
+   also return the radius value for that type
 ------------------------------------------------------------------------- */
 
-void Atom::set_shape(int narg, char **arg)
+int Atom::shape_consistency(int itype,
+			    double &shapex, double &shapey, double &shapez)
 {
-  if (shape == NULL) error->all("Cannot set shape for this atom style");
+  double one[3];
-
+  one[0] = one[1] = one[2] = -1.0;
-  int lo,hi;
+  int flag = 0;
-  force->bounds(arg[0],ntypes,lo,hi);
+  for (int i = 0; i < nlocal; i++) {
-  if (lo < 1 || hi > ntypes) 
+    if (type[i] != itype) continue;
-	error->all("Invalid type for shape set");
+    if (one[0] < 0.0) {
-
+      one[0] = shape[i][0];
-  // store shape as radius, though specified as diameter
+      one[1] = shape[i][1];
-
+      one[2] = shape[i][2];
-  for (int itype = lo; itype <= hi; itype++) {
+    } else if (one[0] != shape[i][0] || one[1] != shape[i][1] || 
-    shape[itype][0] = 0.5*atof(arg[1]);
+	       one[2] != shape[i][2]) flag = 1;
-    shape[itype][1] = 0.5*atof(arg[2]);
-    shape[itype][2] = 0.5*atof(arg[3]);
-    shape_setflag[itype] = 1;
-
-    if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || 
-	shape[itype][2] < 0.0)
-      error->all("Invalid shape value");
-    if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || 
-	shape[itype][2] > 0.0) {
-      if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || 
-	  shape[itype][2] == 0.0)
-	error->all("Invalid shape value");
-    }
   }
-}
 
-/* ----------------------------------------------------------------------
+  int flagall;
-   set all particle shapes as read in from restart file
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
-------------------------------------------------------------------------- */
+  if (flagall) return 0;
 
-void Atom::set_shape(double **values)
+  double oneall[3];
-{
+  MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world);
-  for (int itype = 1; itype <= ntypes; itype++) {
+  shapex = oneall[0];
-    shape[itype][0] = values[itype][0];
+  shapey = oneall[1];
-    shape[itype][1] = values[itype][1];
+  shapez = oneall[2];
-    shape[itype][2] = values[itype][2];
+  return 1;
-    shape_setflag[itype] = 1;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   check that all particle shapes have been set
-------------------------------------------------------------------------- */
-
-void Atom::check_shape()
-{
-  if (shape == NULL) return;
-  for (int itype = 1; itype <= ntypes; itype++)
-    if (shape_setflag[itype] == 0) error->all("All shapes are not set");
-}
-
-/* ----------------------------------------------------------------------
-   set a dipole moment and flag it as set
-   called from reading of data file
-------------------------------------------------------------------------- */
-
-void Atom::set_dipole(const char *str)
-{
-  if (dipole == NULL) error->all("Cannot set dipole for this atom style");
-
-  int itype;
-  double dipole_one;
-  int n = sscanf(str,"%d %lg",&itype,&dipole_one);
-  if (n != 2) error->all("Invalid dipole line in data file");
-
-  dipole[itype] = dipole_one;
-  dipole_setflag[itype] = 1;
-
-  if (dipole[itype] < 0.0) error->all("Invalid dipole value");
-}
-
-/* ----------------------------------------------------------------------
-   set one or more dipole moments and flag them as set
-   called from reading of input script
-------------------------------------------------------------------------- */
-
-void Atom::set_dipole(int narg, char **arg)
-{
-  if (dipole == NULL) error->all("Cannot set dipole for this atom style");
-
-  int lo,hi;
-  force->bounds(arg[0],ntypes,lo,hi);
-  if (lo < 1 || hi > ntypes) error->all("Invalid type for dipole set");
-
-  for (int itype = lo; itype <= hi; itype++) {
-    dipole[itype] = atof(arg[1]);
-    dipole_setflag[itype] = 1;
-
-    if (dipole[itype] < 0.0) error->all("Invalid dipole value");
-  }
-}
-
-/* ----------------------------------------------------------------------
-   set all dipole moments as read in from restart file
-------------------------------------------------------------------------- */
-
-void Atom::set_dipole(double *values)
-{
-  for (int itype = 1; itype <= ntypes; itype++) {
-    dipole[itype] = values[itype];
-    dipole_setflag[itype] = 1;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   check that all dipole moments have been set
-------------------------------------------------------------------------- */
-
-void Atom::check_dipole()
-{
-  if (dipole == NULL) return;
-  for (int itype = 1; itype <= ntypes; itype++)
-    if (dipole_setflag[itype] == 0)
-      error->all("All dipole moments are not set");
 }
 
 /* ----------------------------------------------------------------------

src/atom.h

@@ -48,8 +48,8 @@ class Atom : protected Pointers {
 
   int *molecule;
   double *q,**mu;
-  double **quat,**omega,**angmom,**torque;
+  double **quat,**omega,**angmom,**torque,**shape;
-  double *radius,*density,*rmass,*vfrac,*s0;
+  double *radius,*rmass,*vfrac,*s0;
   double **x0;
 
   int *spin;
@@ -75,15 +75,15 @@ class Atom : protected Pointers {
   int **improper_type;
   int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
 
-  // per-atom array existence flags
+  // atom style and per-atom array existence flags
-  // these can be checked before array is allocated
   // customize by adding new flag
 
-  int molecule_flag;
+  int sphere_flag,ellipsoid_flag,peri_flag,dipole_flag,electron_flag;
-  int q_flag,mu_flag;
+
-  int quat_flag,omega_flag,angmom_flag,torque_flag;
+  int molecule_flag,q_flag,mu_flag;
-  int radius_flag,density_flag,rmass_flag,vfrac_flag;
+  int rmass_flag,radius_flag,omega_flag,torque_flag;
-  int spin_flag,eradius_flag,ervel_flag,erforce_flag;
+  int quat_flag,shape_flag,angmom_flag;
+  int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag;
 
   // extra peratom info in restart file destined for fix & diag 
 
@@ -91,8 +91,8 @@ class Atom : protected Pointers {
 
   // per-type arrays
 
-  double *mass,**shape,*dipole;
+  double *mass;
-  int *mass_setflag,*shape_setflag,*dipole_setflag;
+  int *mass_setflag;
 
   // callback ptrs for atom arrays managed by fix classes
 
@@ -141,14 +141,9 @@ class Atom : protected Pointers {
   void set_mass(int, char **);
   void set_mass(double *);
   void check_mass();
-  void set_shape(const char *);
+
-  void set_shape(int, char **);
+  int radius_consistency(int, double &);
-  void set_shape(double **);
+  int shape_consistency(int, double &, double &, double &);
-  void check_shape();
-  void set_dipole(const char *);
-  void set_dipole(int, char **);
-  void set_dipole(double *);
-  void check_dipole();
 
   void first_reorder();
   void sort();

src/atom_vec.cpp

@@ -24,7 +24,7 @@ AtomVec::AtomVec(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 {
   nmax = 0;
   bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 0;
-  mass_type = shape_type = dipole_type = 0;
+  mass_type = dipole_type = 0;
 }
 
 /* ----------------------------------------------------------------------

src/atom_vec.h

@@ -24,7 +24,6 @@ class AtomVec : protected Pointers {
   int bonds_allow,angles_allow;        // 1 if bonds, angles are used
   int dihedrals_allow,impropers_allow; // 1 if dihedrals, impropers used
   int mass_type;                       // 1 if per-type masses
-  int shape_type;                      // 1 if per-type shape array
   int dipole_type;                     // 1 if per-type dipole moments
 
   int comm_x_only;                     // 1 if only exchange x in forward comm

src/atom_vec_hybrid.cpp

@@ -75,7 +75,6 @@ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp, int narg, char **arg) :
     dihedrals_allow = MAX(dihedrals_allow,styles[k]->dihedrals_allow);
     impropers_allow = MAX(impropers_allow,styles[k]->impropers_allow);
     mass_type = MAX(mass_type,styles[k]->mass_type);
-    shape_type = MAX(shape_type,styles[k]->shape_type);
     dipole_type = MAX(dipole_type,styles[k]->dipole_type);
 
     comm_x_only = MIN(comm_x_only,styles[k]->comm_x_only);

src/atom_vec_sphere.cpp

@@ -15,7 +15,7 @@
 #include "math.h"
 #include "stdlib.h"
 #include "string.h"
-#include "atom_vec_granular.h"
+#include "atom_vec_sphere.h"
 #include "atom.h"
 #include "domain.h"
 #include "modify.h"
@@ -31,7 +31,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) :
+AtomVecSphere::AtomVecSphere(LAMMPS *lmp, int narg, char **arg) :
   AtomVec(lmp, narg, arg)
 {
   molecular = 0;
@@ -46,26 +46,31 @@ AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) :
   size_data_vel = 7;
   xcol_data = 5;
 
-  atom->radius_flag = atom->density_flag = atom->rmass_flag = 1;
+  atom->sphere_flag = 1;
-  atom->omega_flag = atom->torque_flag = 1;
+  atom->radius_flag = atom->rmass_flag = atom->omega_flag = 
+    atom->torque_flag = 1;
 
   PI = 4.0*atan(1.0);
 }
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::init()
+void AtomVecSphere::init()
 {
   AtomVec::init();
 
   // set radvary if particle diameters are time-varying due to fix adapt
 
   radvary = 0;
+  comm_x_only = 1;
+  size_forward = 3;
+
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"adapt") == 0) {
       FixAdapt *fix = (FixAdapt *) modify->fix[i];
       if (fix->diamflag) {
         radvary = 1;
+	comm_x_only = 0;
         size_forward = 5;
       }
     }
@@ -77,7 +82,7 @@ void AtomVecGranular::init()
    n > 0 allocates arrays to size n
 ------------------------------------------------------------------------- */
 
-void AtomVecGranular::grow(int n)
+void AtomVecSphere::grow(int n)
 {
   if (n == 0) nmax += DELTA;
   else nmax = n;
@@ -94,9 +99,7 @@ void AtomVecGranular::grow(int n)
   f = memory->grow(atom->f,nmax,3,"atom:f");
 
   radius = memory->grow(atom->radius,nmax,"atom:radius");
-  density = memory->grow(atom->density,nmax,"atom:density");
   rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
-
   omega = memory->grow(atom->omega,nmax,3,"atom:omega");
   torque = memory->grow(atom->torque,nmax,3,"atom:torque");
 
@@ -109,18 +112,18 @@ void AtomVecGranular::grow(int n)
    reset local array ptrs
 ------------------------------------------------------------------------- */
 
-void AtomVecGranular::grow_reset()
+void AtomVecSphere::grow_reset()
 {
   tag = atom->tag; type = atom->type;
   mask = atom->mask; image = atom->image;
   x = atom->x; v = atom->v; f = atom->f;
-  radius = atom->radius; density = atom->density; rmass = atom->rmass;
+  radius = atom->radius; rmass = atom->rmass;
   omega = atom->omega; torque = atom->torque;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::copy(int i, int j)
+void AtomVecSphere::copy(int i, int j)
 {
   tag[j] = tag[i];
   type[j] = type[i];
@@ -134,7 +137,6 @@ void AtomVecGranular::copy(int i, int j)
   v[j][2] = v[i][2];
 
   radius[j] = radius[i];
-  density[j] = density[i];
   rmass[j] = rmass[i];
   omega[j][0] = omega[i][0];
   omega[j][1] = omega[i][1];
@@ -147,7 +149,7 @@ void AtomVecGranular::copy(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_comm(int n, int *list, double *buf,
+int AtomVecSphere::pack_comm(int n, int *list, double *buf,
 			       int pbc_flag, int *pbc)
 {
   int i,j,m;
@@ -217,7 +219,7 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf,
+int AtomVecSphere::pack_comm_vel(int n, int *list, double *buf,
 				   int pbc_flag, int *pbc)
 {
   int i,j,m;
@@ -361,7 +363,7 @@ int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_comm_one(int i, double *buf)
+int AtomVecSphere::pack_comm_one(int i, double *buf)
 {
   if (radvary == 0) return 0;
 
@@ -372,7 +374,7 @@ int AtomVecGranular::pack_comm_one(int i, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::unpack_comm(int n, int first, double *buf)
+void AtomVecSphere::unpack_comm(int n, int first, double *buf)
 {
   int i,m,last;
 
@@ -399,7 +401,7 @@ void AtomVecGranular::unpack_comm(int n, int first, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf)
+void AtomVecSphere::unpack_comm_vel(int n, int first, double *buf)
 {
   int i,m,last;
 
@@ -438,7 +440,7 @@ void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::unpack_comm_one(int i, double *buf)
+int AtomVecSphere::unpack_comm_one(int i, double *buf)
 {
   if (radvary == 0) return 0;
 
@@ -449,7 +451,7 @@ int AtomVecGranular::unpack_comm_one(int i, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_reverse(int n, int first, double *buf)
+int AtomVecSphere::pack_reverse(int n, int first, double *buf)
 {
   int i,m,last;
 
@@ -468,7 +470,7 @@ int AtomVecGranular::pack_reverse(int n, int first, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_reverse_one(int i, double *buf)
+int AtomVecSphere::pack_reverse_one(int i, double *buf)
 {
   buf[0] = torque[i][0];
   buf[1] = torque[i][1];
@@ -478,7 +480,7 @@ int AtomVecGranular::pack_reverse_one(int i, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::unpack_reverse(int n, int *list, double *buf)
+void AtomVecSphere::unpack_reverse(int n, int *list, double *buf)
 {
   int i,j,m;
 
@@ -496,7 +498,7 @@ void AtomVecGranular::unpack_reverse(int n, int *list, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::unpack_reverse_one(int i, double *buf)
+int AtomVecSphere::unpack_reverse_one(int i, double *buf)
 {
   torque[i][0] += buf[0];
   torque[i][1] += buf[1];
@@ -506,7 +508,7 @@ int AtomVecGranular::unpack_reverse_one(int i, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_border(int n, int *list, double *buf,
+int AtomVecSphere::pack_border(int n, int *list, double *buf,
 				 int pbc_flag, int *pbc)
 {
   int i,j,m;
@@ -552,7 +554,7 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_border_vel(int n, int *list, double *buf,
+int AtomVecSphere::pack_border_vel(int n, int *list, double *buf,
 				     int pbc_flag, int *pbc)
 {
   int i,j,m;
@@ -639,7 +641,7 @@ int AtomVecGranular::pack_border_vel(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_border_one(int i, double *buf)
+int AtomVecSphere::pack_border_one(int i, double *buf)
 {
   buf[0] = radius[i];
   buf[1] = rmass[i];
@@ -648,7 +650,7 @@ int AtomVecGranular::pack_border_one(int i, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::unpack_border(int n, int first, double *buf)
+void AtomVecSphere::unpack_border(int n, int first, double *buf)
 {
   int i,m,last;
 
@@ -670,7 +672,7 @@ void AtomVecGranular::unpack_border(int n, int first, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-void AtomVecGranular::unpack_border_vel(int n, int first, double *buf)
+void AtomVecSphere::unpack_border_vel(int n, int first, double *buf)
 {
   int i,m,last;
 
@@ -697,7 +699,7 @@ void AtomVecGranular::unpack_border_vel(int n, int first, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::unpack_border_one(int i, double *buf)
+int AtomVecSphere::unpack_border_one(int i, double *buf)
 {
   radius[i] = buf[0];
   rmass[i] = buf[1];
@@ -709,7 +711,7 @@ int AtomVecGranular::unpack_border_one(int i, double *buf)
    xyz must be 1st 3 values, so comm::exchange() can test on them 
 ------------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_exchange(int i, double *buf)
+int AtomVecSphere::pack_exchange(int i, double *buf)
 {
   int m = 1;
   buf[m++] = x[i][0];
@@ -724,7 +726,6 @@ int AtomVecGranular::pack_exchange(int i, double *buf)
   buf[m++] = image[i];
 
   buf[m++] = radius[i];
-  buf[m++] = density[i];
   buf[m++] = rmass[i];
   buf[m++] = omega[i][0];
   buf[m++] = omega[i][1];
@@ -740,7 +741,7 @@ int AtomVecGranular::pack_exchange(int i, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecGranular::unpack_exchange(double *buf)
+int AtomVecSphere::unpack_exchange(double *buf)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -758,7 +759,6 @@ int AtomVecGranular::unpack_exchange(double *buf)
   image[nlocal] = static_cast<int> (buf[m++]);
 
   radius[nlocal] = buf[m++];
-  density[nlocal] = buf[m++];
   rmass[nlocal] = buf[m++];
   omega[nlocal][0] = buf[m++];
   omega[nlocal][1] = buf[m++];
@@ -778,7 +778,7 @@ int AtomVecGranular::unpack_exchange(double *buf)
    include extra data stored by fixes
 ------------------------------------------------------------------------- */
 
-int AtomVecGranular::size_restart()
+int AtomVecSphere::size_restart()
 {
   int i;
 
@@ -799,7 +799,7 @@ int AtomVecGranular::size_restart()
    molecular types may be negative, but write as positive   
 ------------------------------------------------------------------------- */
 
-int AtomVecGranular::pack_restart(int i, double *buf)
+int AtomVecSphere::pack_restart(int i, double *buf)
 {
   int m = 1;
   buf[m++] = x[i][0];
@@ -814,7 +814,7 @@ int AtomVecGranular::pack_restart(int i, double *buf)
   buf[m++] = v[i][2];
 
   buf[m++] = radius[i];
-  buf[m++] = density[i];
+  buf[m++] = rmass[i];
   buf[m++] = omega[i][0];
   buf[m++] = omega[i][1];
   buf[m++] = omega[i][2];
@@ -831,7 +831,7 @@ int AtomVecGranular::pack_restart(int i, double *buf)
    unpack data for one atom from restart file including extra quantities
 ------------------------------------------------------------------------- */
 
-int AtomVecGranular::unpack_restart(double *buf)
+int AtomVecSphere::unpack_restart(double *buf)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) {
@@ -853,13 +853,7 @@ int AtomVecGranular::unpack_restart(double *buf)
   v[nlocal][2] = buf[m++];
 
   radius[nlocal] = buf[m++];
-  density[nlocal] = buf[m++];
+  rmass[nlocal] = buf[m++];
-
-  if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal];
-  else 
-    rmass[nlocal] = 4.0*PI/3.0 * 
-      radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
-
   omega[nlocal][0] = buf[m++];
   omega[nlocal][1] = buf[m++];
   omega[nlocal][2] = buf[m++];
@@ -879,7 +873,7 @@ int AtomVecGranular::unpack_restart(double *buf)
    set other values to defaults
 ------------------------------------------------------------------------- */
 
-void AtomVecGranular::create_atom(int itype, double *coord)
+void AtomVecSphere::create_atom(int itype, double *coord)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -896,9 +890,7 @@ void AtomVecGranular::create_atom(int itype, double *coord)
   v[nlocal][2] = 0.0;
 
   radius[nlocal] = 0.5;
-  density[nlocal] = 1.0;
+  rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal];
-  rmass[nlocal] = 4.0*PI/3.0 *
-    radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
   omega[nlocal][0] = 0.0;
   omega[nlocal][1] = 0.0;
   omega[nlocal][2] = 0.0;
@@ -911,7 +903,7 @@ void AtomVecGranular::create_atom(int itype, double *coord)
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values)
+void AtomVecSphere::data_atom(double *coord, int imagetmp, char **values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -928,14 +920,14 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values)
   if (radius[nlocal] < 0.0)
     error->one("Invalid radius in Atoms section of data file");
 
-  density[nlocal] = atof(values[3]);
+  double density = atof(values[3]);
-  if (density[nlocal] <= 0.0)
+  if (density <= 0.0)
     error->one("Invalid density in Atoms section of data file");
 
-  if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal];
+  if (radius[nlocal] == 0.0) rmass[nlocal] = density;
   else 
     rmass[nlocal] = 4.0*PI/3.0 *
-      radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
+      radius[nlocal]*radius[nlocal]*radius[nlocal] * density;
 
   x[nlocal][0] = coord[0];
   x[nlocal][1] = coord[1];
@@ -959,20 +951,20 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values)
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecGranular::data_atom_hybrid(int nlocal, char **values)
+int AtomVecSphere::data_atom_hybrid(int nlocal, char **values)
 {
   radius[nlocal] = 0.5 * atof(values[0]);
   if (radius[nlocal] < 0.0)
     error->one("Invalid radius in Atoms section of data file");
 
-  density[nlocal] = atof(values[1]);
+  double density = atof(values[1]);
-  if (density[nlocal] <= 0.0)
+  if (density <= 0.0)
     error->one("Invalid density in Atoms section of data file");
 
-  if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal];
+  if (radius[nlocal] == 0.0) rmass[nlocal] = density;
   else 
     rmass[nlocal] = 4.0*PI/3.0 *
-      radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
+      radius[nlocal]*radius[nlocal]*radius[nlocal] * density;
 
   return 2;
 }
@@ -981,7 +973,7 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values)
    unpack one line from Velocities section of data file
 ------------------------------------------------------------------------- */
 
-void AtomVecGranular::data_vel(int m, char **values)
+void AtomVecSphere::data_vel(int m, char **values)
 {
   v[m][0] = atof(values[0]);
   v[m][1] = atof(values[1]);
@@ -995,7 +987,7 @@ void AtomVecGranular::data_vel(int m, char **values)
    unpack hybrid quantities from one line in Velocities section of data file
 ------------------------------------------------------------------------- */
 
-int AtomVecGranular::data_vel_hybrid(int m, char **values)
+int AtomVecSphere::data_vel_hybrid(int m, char **values)
 {
   omega[m][0] = atof(values[0]);
   omega[m][1] = atof(values[1]);
@@ -1007,7 +999,7 @@ int AtomVecGranular::data_vel_hybrid(int m, char **values)
    return # of bytes of allocated memory 
 ------------------------------------------------------------------------- */
 
-bigint AtomVecGranular::memory_usage()
+bigint AtomVecSphere::memory_usage()
 {
   bigint bytes = 0;
 
@@ -1020,7 +1012,6 @@ bigint AtomVecGranular::memory_usage()
   if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
 
   if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax);
-  if (atom->memcheck("density")) bytes += memory->usage(density,nmax);
   if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
   if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
   if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);

src/atom_vec_sphere.h

@@ -13,21 +13,21 @@
 
 #ifdef ATOM_CLASS
 
-AtomStyle(granular,AtomVecGranular)
+AtomStyle(sphere,AtomVecSphere)
 
 #else
 
-#ifndef LMP_ATOM_VEC_GRANULAR_H
+#ifndef LMP_ATOM_VEC_SPHERE_H
-#define LMP_ATOM_VEC_GRANULAR_H
+#define LMP_ATOM_VEC_SPHERE_H
 
 #include "atom_vec.h"
 
 namespace LAMMPS_NS {
 
-class AtomVecGranular : public AtomVec {
+class AtomVecSphere : public AtomVec {
  public:
-  AtomVecGranular(class LAMMPS *, int, char **);
+  AtomVecSphere(class LAMMPS *, int, char **);
-  ~AtomVecGranular() {}
+  ~AtomVecSphere() {}
   void init();
   void grow(int);
   void grow_reset();

src/compute_erotate_sphere.cpp

@@ -35,40 +35,16 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   extscalar = 1;
 
-  // error checks
+  // error check
 
-  if (!atom->omega_flag) 
+  if (!atom->sphere_flag) 
-    error->all("Compute erotate/sphere requires atom attribute omega");
+    error->all("Compute erotate/sphere requires atom style sphere");
-  if (!atom->radius_flag && !atom->avec->shape_type)
-    error->all("Compute erotate/sphere requires atom attribute "
-	       "radius or shape");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeERotateSphere::init()
 {
-  int i,itype;
-
-  // if shape used, check that all particles are spherical
-  // point particles are allowed
-
-  if (atom->radius == NULL) {
-    double **shape = atom->shape;
-    int *type = atom->type;
-    int *mask = atom->mask;
-    int nlocal = atom->nlocal;
-
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	itype = type[i];
-	if (shape[itype][0] != shape[itype][1] || 
-	    shape[itype][0] != shape[itype][2])
-	  error->one("Compute erotate/sphere requires "
-		     "spherical particle shapes");
-      }
-  }
-
   pfactor = 0.5 * force->mvv2e * INERTIA;
 }
 
@@ -76,59 +52,22 @@ void ComputeERotateSphere::init()
 
 double ComputeERotateSphere::compute_scalar()
 {
-  int i,itype;
-
   invoked_scalar = update->ntimestep;
 
   double **omega = atom->omega;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  double **shape = atom->shape;
   int *mask = atom->mask;
-  int *type = atom->type;
   int nlocal = atom->nlocal;
 
   // sum rotational energy for each particle
-  // point particles will not contribute due to radius or shape = 0
+  // point particles will not contribute due to radius = 0
 
   double erotate = 0.0;
-
+  for (int i = 0; i < nlocal; i++) 
-  if (radius) {
+    if (mask[i] & groupbit)
-    if (rmass) {
+      erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-      for (i = 0; i < nlocal; i++) 
+		  omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];
-	if (mask[i] & groupbit)
-	  erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		      omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];
-    } else {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		      omega[i][2]*omega[i][2]) * 
-	    radius[i]*radius[i]*mass[itype];
-	}
-    }
-    
-  } else {
-    if (rmass) {
-      for (i = 0; i < nlocal; i++) 
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		      omega[i][2]*omega[i][2]) * 
-	    shape[itype][0]*shape[itype][0]*rmass[i];
-	}
-    } else {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		      omega[i][2]*omega[i][2]) *
-	    shape[itype][0]*shape[itype][0]*mass[itype];
-	}
-    }
-  }
 
   MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   scalar *= pfactor;

src/compute_property_atom.cpp

@@ -124,6 +124,11 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
 	error->all("Compute property/atom for "
 		   "atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_muz;
+    } else if (strcmp(arg[iarg],"mu") == 0) {
+      if (!atom->mu_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_mu;
 
     } else if (strcmp(arg[iarg],"radius") == 0) {
       if (!atom->radius_flag)
@@ -161,6 +166,21 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
 		   "atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_angmomz;
 
+    } else if (strcmp(arg[iarg],"shapex") == 0) {
+      if (!atom->shape_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_shapex;
+    } else if (strcmp(arg[iarg],"shapey") == 0) {
+      if (!atom->shape_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_shapey;
+    } else if (strcmp(arg[iarg],"shapez") == 0) {
+      if (!atom->shape_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_shapez;
     } else if (strcmp(arg[iarg],"quatw") == 0) {
       if (!atom->quat_flag)
 	error->all("Compute property/atom for "
@@ -831,6 +851,21 @@ void ComputePropertyAtom::pack_muz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void ComputePropertyAtom::pack_mu(int n)
+{
+  double **mu = atom->mu;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = mu[i][3];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ComputePropertyAtom::pack_radius(int n)
 {
   double *radius = atom->radius;
@@ -936,6 +971,51 @@ void ComputePropertyAtom::pack_angmomz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void ComputePropertyAtom::pack_shapex(int n)
+{
+  double **shape = atom->shape;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = shape[i][0];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_shapey(int n)
+{
+  double **shape = atom->shape;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = shape[i][1];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_shapez(int n)
+{
+  double **shape = atom->shape;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = shape[i][2];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ComputePropertyAtom::pack_quatw(int n)
 {
   double **quat = atom->quat;

src/compute_property_atom.h

@@ -76,6 +76,7 @@ class ComputePropertyAtom : public Compute {
   void pack_mux(int);
   void pack_muy(int);
   void pack_muz(int);
+  void pack_mu(int);
   void pack_radius(int);
   void pack_omegax(int);
   void pack_omegay(int);
@@ -83,6 +84,9 @@ class ComputePropertyAtom : public Compute {
   void pack_angmomx(int);
   void pack_angmomy(int);
   void pack_angmomz(int);
+  void pack_shapex(int);
+  void pack_shapey(int);
+  void pack_shapez(int);
   void pack_quatw(int);
   void pack_quati(int);
   void pack_quatj(int);

src/compute_temp_sphere.cpp

@@ -55,11 +55,8 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
 
   // error checks
 
-  if (!atom->omega_flag)
+  if (!atom->sphere_flag) 
-    error->all("Compute temp/sphere requires atom attribute omega");
+    error->all("Compute temp/sphere requires atom style sphere");
-  if (!atom->radius_flag && !atom->avec->shape_type)
-    error->all("Compute temp/sphere requires atom attribute "
-	       "radius or shape");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -74,29 +71,8 @@ ComputeTempSphere::~ComputeTempSphere()
 
 void ComputeTempSphere::init()
 {
-  int i,itype;
-
-  // if shape used, check that all particles are spherical
-  // point particles are allowed
-
-  if (atom->radius == NULL) {
-    double **shape = atom->shape;
-    int *type = atom->type;
-    int *mask = atom->mask;
-    int nlocal = atom->nlocal;
-
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	itype = type[i];
-	if (shape[itype][0] != shape[itype][1] || 
-	    shape[itype][0] != shape[itype][2])
-	  error->one("Compute temp/sphere requires "
-		     "spherical particle shapes");
-      }
-  }
-
   if (tempbias) {
-    i = modify->find_compute(id_bias);
+    int i = modify->find_compute(id_bias);
     if (i < 0) error->all("Could not find compute ID for temperature bias");
     tbias = modify->compute[i];
     if (tbias->tempflag == 0)
@@ -111,7 +87,7 @@ void ComputeTempSphere::init()
   }
 
   fix_dof = 0;
-  for (i = 0; i < modify->nfix; i++)
+  for (int i = 0; i < modify->nfix; i++)
     fix_dof += modify->fix[i]->dof(igroup);
   dof_compute();
 }
@@ -130,44 +106,22 @@ void ComputeTempSphere::dof_compute()
   int dimension = domain->dimension;
 
   double *radius = atom->radius;
-  double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   count = 0;
   if (dimension == 3) {
-    if (radius) {
+    for (int i = 0; i < nlocal; i++)
-      for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
-	if (mask[i] & groupbit) {
+	if (radius[i] == 0.0) count += 3;
-	  if (radius[i] == 0.0) count += 3;
+	else count += 6;
-	  else count += 6;
-	}
       }
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  if (shape[type[i]][0] == 0.0) count += 3;
-	  else count += 6;
-	}
-      }
-    }
   } else {
-    if (radius) {
+    for (int i = 0; i < nlocal; i++)
-      for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
-	if (mask[i] & groupbit) {
+	if (radius[i] == 0.0) count += 2;
-	  if (radius[i] == 0.0) count += 2;
+	else count += 3;
-	  else count += 3;
-	}
       }
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  if (shape[type[i]][0] == 0.0) count += 2;
-	  else count += 3;
-	}
-      }
-    }
   }
 
   MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world);
@@ -185,41 +139,21 @@ void ComputeTempSphere::dof_compute()
 
     count = 0;
     if (dimension == 3) {
-      if (radius) {
+      for (int i = 0; i < nlocal; i++)
-	for (int i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) {
-	  if (mask[i] & groupbit) {
+	  if (tbias->dof_remove(i)) {
-	    if (tbias->dof_remove(i)) {
+	    if (radius[i] == 0.0) count += 3;
-	      if (radius[i] == 0.0) count += 3;
+	    else count += 6;
-	      else count += 6;
-	    }
 	  }
-      } else {
+	}
-	for (int i = 0; i < nlocal; i++)
-	  if (mask[i] & groupbit) {
-	    if (tbias->dof_remove(i)) {
-	      if (shape[type[i]][0] == 0.0) count += 3;
-	      else count += 6;
-	    }
-	  }
-      }
     } else {
-      if (radius) {
+      for (int i = 0; i < nlocal; i++)
-	for (int i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) {
-	  if (mask[i] & groupbit) {
+	  if (tbias->dof_remove(i)) {
-	    if (tbias->dof_remove(i)) {
+	    if (radius[i] == 0.0) count += 2;
-	      if (radius[i] == 0.0) count += 2;
+	    else count += 3;
-	      else count += 3;
-	    }
 	  }
-      } else {
+	}
-	for (int i = 0; i < nlocal; i++)
-	  if (mask[i] & groupbit) {
-	    if (tbias->dof_remove(i)) {
-	      if (shape[type[i]][0] == 0.0) count += 2;
-	      else count += 3;
-	    }
-	  }
-      }
     }
 
     MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world);
@@ -235,8 +169,6 @@ void ComputeTempSphere::dof_compute()
 
 double ComputeTempSphere::compute_scalar()
 {
-  int i,itype;
-
   invoked_scalar = update->ntimestep;
 
   if (tempbias) {
@@ -248,65 +180,20 @@ double ComputeTempSphere::compute_scalar()
   double **omega = atom->omega;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  // 4 cases depending on radius vs shape and rmass vs mass
+  // point particles will not contribute rotation due to radius = 0
-  // point particles will not contribute rotation due to radius or shape = 0
 
   double t = 0.0;
 
-  if (radius) {
+  for (int i = 0; i < nlocal; i++)
-    if (rmass) {
+    if (mask[i] & groupbit) {
-      for (i = 0; i < nlocal; i++)
+      t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
-	if (mask[i] & groupbit) {
+      t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-	  t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * 
+	    omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i];
-	    rmass[i];
-	  t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		omega[i][2]*omega[i][2]) * 
-	    INERTIA*radius[i]*radius[i]*rmass[i];
-	}
-
-    } else {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * 
-	    mass[itype];
-	  t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		omega[i][2]*omega[i][2]) * 
-	    INERTIA*radius[i]*radius[i]*mass[itype];
-	}
     }
 
-  } else {
-    if (rmass) {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * 
-	    rmass[i];
-	  t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		omega[i][2]*omega[i][2]) *
-	    INERTIA*shape[itype][0]*shape[itype][0]*rmass[i];
-	}
-      
-    } else {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * 
-	    mass[itype];
-	  t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-		omega[i][2]*omega[i][2]) *
-	    INERTIA*shape[itype][0]*shape[itype][0]*mass[itype];
-	}
-    }
-  }
-
   if (tempbias) tbias->restore_bias_all();
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
@@ -319,8 +206,6 @@ double ComputeTempSphere::compute_scalar()
 
 void ComputeTempSphere::compute_vector()
 {
-  int i,itype;
-
   invoked_vector = update->ntimestep;
 
   if (tempbias) {
@@ -330,112 +215,39 @@ void ComputeTempSphere::compute_vector()
 
   double **v = atom->v;
   double **omega = atom->omega;
-  double *mass = atom->mass;
   double *rmass = atom->rmass;
   double *radius = atom->radius;
-  double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  // 4 cases depending on radius vs shape and rmass vs mass
+  // point particles will not contribute rotation due to radius = 0
-  // point particles will not contribute rotation due to radius or shape = 0
 
   double massone,inertiaone,t[6];
-  for (i = 0; i < 6; i++) t[i] = 0.0;
+  for (int i = 0; i < 6; i++) t[i] = 0.0;
 
-  if (radius) {
+  for (int i = 0; i < nlocal; i++)
-    if (rmass) {
+    if (mask[i] & groupbit) {
-      for (i = 0; i < nlocal; i++)
+      massone = rmass[i];
-	if (mask[i] & groupbit) {
+      t[0] += massone * v[i][0]*v[i][0];
-	  massone = rmass[i];
+      t[1] += massone * v[i][1]*v[i][1];
-	  t[0] += massone * v[i][0]*v[i][0];
+      t[2] += massone * v[i][2]*v[i][2];
-	  t[1] += massone * v[i][1]*v[i][1];
+      t[3] += massone * v[i][0]*v[i][1];
-	  t[2] += massone * v[i][2]*v[i][2];
+      t[4] += massone * v[i][0]*v[i][2];
-	  t[3] += massone * v[i][0]*v[i][1];
+      t[5] += massone * v[i][1]*v[i][2];
-	  t[4] += massone * v[i][0]*v[i][2];
+      
-	  t[5] += massone * v[i][1]*v[i][2];
+      inertiaone = INERTIA*radius[i]*radius[i]*rmass[i];
-	  
+      t[0] += inertiaone * omega[i][0]*omega[i][0];
-	  inertiaone = INERTIA*radius[i]*radius[i]*rmass[i];
+      t[1] += inertiaone * omega[i][1]*omega[i][1];
-	  t[0] += inertiaone * omega[i][0]*omega[i][0];
+      t[2] += inertiaone * omega[i][2]*omega[i][2];
-	  t[1] += inertiaone * omega[i][1]*omega[i][1];
+      t[3] += inertiaone * omega[i][0]*omega[i][1];
-	  t[2] += inertiaone * omega[i][2]*omega[i][2];
+      t[4] += inertiaone * omega[i][0]*omega[i][2];
-	  t[3] += inertiaone * omega[i][0]*omega[i][1];
+      t[5] += inertiaone * omega[i][1]*omega[i][2];
-	  t[4] += inertiaone * omega[i][0]*omega[i][2];
-	  t[5] += inertiaone * omega[i][1]*omega[i][2];
-	}
-
-    } else {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  massone = mass[itype];
-	  t[0] += massone * v[i][0]*v[i][0];
-	  t[1] += massone * v[i][1]*v[i][1];
-	  t[2] += massone * v[i][2]*v[i][2];
-	  t[3] += massone * v[i][0]*v[i][1];
-	  t[4] += massone * v[i][0]*v[i][2];
-	  t[5] += massone * v[i][1]*v[i][2];
-	  
-	  inertiaone = INERTIA*radius[i]*radius[i]*mass[itype];
-	  t[0] += inertiaone * omega[i][0]*omega[i][0];
-	  t[1] += inertiaone * omega[i][1]*omega[i][1];
-	  t[2] += inertiaone * omega[i][2]*omega[i][2];
-	  t[3] += inertiaone * omega[i][0]*omega[i][1];
-	  t[4] += inertiaone * omega[i][0]*omega[i][2];
-	  t[5] += inertiaone * omega[i][1]*omega[i][2];
-	}
     }
 
-  } else {
-    if (rmass) {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  massone = rmass[i];
-	  t[0] += massone * v[i][0]*v[i][0];
-	  t[1] += massone * v[i][1]*v[i][1];
-	  t[2] += massone * v[i][2]*v[i][2];
-	  t[3] += massone * v[i][0]*v[i][1];
-	  t[4] += massone * v[i][0]*v[i][2];
-	  t[5] += massone * v[i][1]*v[i][2];
-	  
-	  inertiaone = INERTIA*shape[itype][0]*shape[itype][0]*rmass[i];
-	  t[0] += inertiaone * omega[i][0]*omega[i][0];
-	  t[1] += inertiaone * omega[i][1]*omega[i][1];
-	  t[2] += inertiaone * omega[i][2]*omega[i][2];
-	  t[3] += inertiaone * omega[i][0]*omega[i][1];
-	  t[4] += inertiaone * omega[i][0]*omega[i][2];
-	  t[5] += inertiaone * omega[i][1]*omega[i][2];
-	}
-
-    } else {
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  massone = mass[itype];
-	  t[0] += massone * v[i][0]*v[i][0];
-	  t[1] += massone * v[i][1]*v[i][1];
-	  t[2] += massone * v[i][2]*v[i][2];
-	  t[3] += massone * v[i][0]*v[i][1];
-	  t[4] += massone * v[i][0]*v[i][2];
-	  t[5] += massone * v[i][1]*v[i][2];
-	  
-	  inertiaone = INERTIA*shape[itype][0]*shape[itype][0]*mass[itype];
-	  t[0] += inertiaone * omega[i][0]*omega[i][0];
-	  t[1] += inertiaone * omega[i][1]*omega[i][1];
-	  t[2] += inertiaone * omega[i][2]*omega[i][2];
-	  t[3] += inertiaone * omega[i][0]*omega[i][1];
-	  t[4] += inertiaone * omega[i][0]*omega[i][2];
-	  t[5] += inertiaone * omega[i][1]*omega[i][2];
-	}
-    }
-  }
-
   if (tempbias) tbias->restore_bias_all();
 
   MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
+  for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
 }
 
 /* ----------------------------------------------------------------------

src/dump_custom.cpp

@@ -35,11 +35,12 @@ using namespace LAMMPS_NS;
 // same list as in compute_property.cpp, also customize that command
 
 enum{ID,MOL,TYPE,MASS,
-       X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ,
+     X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ,
-       VX,VY,VZ,FX,FY,FZ,
+     VX,VY,VZ,FX,FY,FZ,
-       Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
+     Q,MUX,MUY,MUZ,MU,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
-       QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
+     SHAPEX,SHAPEY,SHAPEZ,
-       COMPUTE,FIX,VARIABLE};
+     QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
+     COMPUTE,FIX,VARIABLE};
 enum{LT,LE,GT,GE,EQ,NEQ};
 enum{INT,DOUBLE};
 
@@ -610,17 +611,22 @@ int DumpCustom::count()
 	if (!atom->mu_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
 	ptr = &atom->mu[0][0];
-	nstride = 3;
+	nstride = 4;
       } else if (thresh_array[ithresh] == MUY) {
 	if (!atom->mu_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
 	ptr = &atom->mu[0][1];
-	nstride = 3;
+	nstride = 4;
       } else if (thresh_array[ithresh] == MUZ) {
 	if (!atom->mu_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
 	ptr = &atom->mu[0][2];
-	nstride = 3;
+	nstride = 4;
+      } else if (thresh_array[ithresh] == MU) {
+	if (!atom->mu_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+	ptr = &atom->mu[0][3];
+	nstride = 4;
 
       } else if (thresh_array[ithresh] == RADIUS) {
 	if (!atom->radius_flag)
@@ -658,6 +664,21 @@ int DumpCustom::count()
 	ptr = &atom->angmom[0][2];
 	nstride = 3;
 
+      } else if (thresh_array[ithresh] == SHAPEX) {
+	if (!atom->shape_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+	ptr = &atom->shape[0][0];
+	nstride = 3;
+      } else if (thresh_array[ithresh] == SHAPEY) {
+	if (!atom->shape_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+	ptr = &atom->shape[0][1];
+	nstride = 3;
+      } else if (thresh_array[ithresh] == SHAPEZ) {
+	if (!atom->shape_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+	ptr = &atom->shape[0][2];
+	nstride = 3;
       } else if (thresh_array[ithresh] == QUATW) {
 	if (!atom->quat_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
@@ -944,6 +965,11 @@ void DumpCustom::parse_fields(int narg, char **arg)
 	error->all("Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_muz;
       vtype[i] = DOUBLE;
+    } else if (strcmp(arg[iarg],"mu") == 0) {
+      if (!atom->mu_flag)
+	error->all("Dumping an atom property that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_mu;
+      vtype[i] = DOUBLE;
 
     } else if (strcmp(arg[iarg],"radius") == 0) {
       if (!atom->radius_flag)
@@ -981,6 +1007,21 @@ void DumpCustom::parse_fields(int narg, char **arg)
       pack_choice[i] = &DumpCustom::pack_angmomz;
       vtype[i] = DOUBLE;
 
+    } else if (strcmp(arg[iarg],"shapex") == 0) {
+      if (!atom->shape_flag)
+	error->all("Dumping an atom property that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_shapex;
+      vtype[i] = DOUBLE;
+    } else if (strcmp(arg[iarg],"shapey") == 0) {
+      if (!atom->shape_flag)
+	error->all("Dumping an atom property that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_shapey;
+      vtype[i] = DOUBLE;
+    } else if (strcmp(arg[iarg],"shapez") == 0) {
+      if (!atom->shape_flag)
+	error->all("Dumping an atom property that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_shapez;
+      vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"quatw") == 0) {
       if (!atom->quat_flag)
 	error->all("Dumping an atom property that isn't allocated");
@@ -1300,6 +1341,7 @@ int DumpCustom::modify_param(int narg, char **arg)
     else if (strcmp(arg[1],"mux") == 0) thresh_array[nthresh] = MUX;
     else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY;
     else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ;
+    else if (strcmp(arg[1],"mu") == 0) thresh_array[nthresh] = MU;
 
     else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
     else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX;
@@ -1309,6 +1351,9 @@ int DumpCustom::modify_param(int narg, char **arg)
     else if (strcmp(arg[1],"angmomy") == 0) thresh_array[nthresh] = ANGMOMY;
     else if (strcmp(arg[1],"angmomz") == 0) thresh_array[nthresh] = ANGMOMZ;
 
+    else if (strcmp(arg[1],"shapex") == 0) thresh_array[nthresh] = SHAPEX;
+    else if (strcmp(arg[1],"shapey") == 0) thresh_array[nthresh] = SHAPEY;
+    else if (strcmp(arg[1],"shapez") == 0) thresh_array[nthresh] = SHAPEZ;
     else if (strcmp(arg[1],"quatw") == 0) thresh_array[nthresh] = QUATW;
     else if (strcmp(arg[1],"quati") == 0) thresh_array[nthresh] = QUATI;
     else if (strcmp(arg[1],"quatj") == 0) thresh_array[nthresh] = QUATJ;
@@ -2039,6 +2084,20 @@ void DumpCustom::pack_muz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::pack_mu(int n)
+{
+  double **mu = atom->mu;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = mu[i][3];
+      n += size_one;
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::pack_radius(int n)
 {
   double *radius = atom->radius;
@@ -2137,6 +2196,48 @@ void DumpCustom::pack_angmomz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::pack_shapex(int n)
+{
+  double **shape = atom->shape;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = shape[i][0];
+      n += size_one;
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_shapey(int n)
+{
+  double **shape = atom->shape;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = shape[i][1];
+      n += size_one;
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_shapez(int n)
+{
+  double **shape = atom->shape;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = shape[i][2];
+      n += size_one;
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::pack_quatw(int n)
 {
   double **quat = atom->quat;

src/dump_custom.h

@@ -134,6 +134,7 @@ class DumpCustom : public Dump {
   void pack_mux(int);
   void pack_muy(int);
   void pack_muz(int);
+  void pack_mu(int);
   void pack_radius(int);
   void pack_omegax(int);
   void pack_omegay(int);
@@ -141,6 +142,9 @@ class DumpCustom : public Dump {
   void pack_angmomx(int);
   void pack_angmomy(int);
   void pack_angmomz(int);
+  void pack_shapex(int);
+  void pack_shapey(int);
+  void pack_shapez(int);
   void pack_quatw(int);
   void pack_quati(int);
   void pack_quatj(int);

src/fix_adapt.cpp

@@ -307,25 +307,31 @@ void FixAdapt::change_settings()
     } else if (ad->which == ATOM) {
 
       // set radius from diameter
-      // also set rmass if both rmass and density are defined
+      // also scale rmass to new value
 
       if (ad->aparam == DIAMETER) {
 	int mflag = 0;
-	if (atom->rmass_flag && atom->density_flag) mflag = 1;
+	if (atom->rmass_flag) mflag = 1;
 	double PI = 4.0*atan(1.0);
+	double density;
 
 	double *radius = atom->radius;
 	double *rmass = atom->rmass;
-	double *density = atom->density;
 	int *mask = atom->mask;
 	int nlocal = atom->nlocal;
 
-	for (i = 0; i < nlocal; i++)
+	if (mflag == 0) {
-	  if (mask[i] & groupbit) {
+	  for (i = 0; i < nlocal; i++)
-	    radius[i] = 0.5*value;
+	    if (mask[i] & groupbit)
-	    if (mflag) rmass[i] = 4.0*PI/3.0 * 
+	      radius[i] = 0.5*value;
-			 radius[i]*radius[i]*radius[i] * density[i];
+	} else {
-	  }
+	  for (i = 0; i < nlocal; i++)
+	    if (mask[i] & groupbit) {
+	      density = rmass[i] / (4.0*PI/3.0 * radius[i]*radius[i]*radius[i]);
+	      radius[i] = 0.5*value;
+	      rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density;
+	    }
+	}
       }
     }
   }

src/fix_nh_sphere.cpp

@@ -31,12 +31,8 @@ using namespace LAMMPS_NS;
 FixNHSphere::FixNHSphere(LAMMPS *lmp, int narg, char **arg) :
   FixNH(lmp, narg, arg)
 {
-  if (!atom->omega_flag || !atom->torque_flag)
+  if (!atom->sphere_flag)
-    error->all("Fix nvt/nph/npt sphere requires "
+    error->all("Fix nvt/nph/npt sphere requires atom style sphere");
-	       "atom attributes omega, torque");
-  if (!atom->radius_flag && !atom->avec->shape_type)
-    error->all("Fix nvt/nph/npt sphere requires "
-	       "atom attribute diameter or shape");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -45,36 +41,17 @@ void FixNHSphere::init()
 {
   int i,itype;
 
-  // check that all particles are finite-size and spherical
+  // check that all particles are finite-size
   // no point particles allowed
 
-  if (atom->radius_flag) {
+  double *radius = atom->radius;
-    double *radius = atom->radius;
+  int *mask = atom->mask;
-    int *mask = atom->mask;
+  int nlocal = atom->nlocal;
-    int nlocal = atom->nlocal;
 
-    for (i = 0; i < nlocal; i++)
+  for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
+    if (mask[i] & groupbit)
-	if (radius[i] == 0.0)
+      if (radius[i] == 0.0)
-	  error->one("Fix nvt/sphere requires extended particles");
+	error->one("Fix nvt/sphere requires extended particles");
-      }
-
-  } else {
-    double **shape = atom->shape;
-    int *type = atom->type;
-    int *mask = atom->mask;
-    int nlocal = atom->nlocal;
-
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	itype = type[i];
-	if (shape[itype][0] == 0.0)
-	  error->one("Fix nvt/sphere requires extended particles");
-	if (shape[itype][0] != shape[itype][1] || 
-	    shape[itype][0] != shape[itype][2])
-	  error->one("Fix nvt/sphere requires spherical particle shapes");
-      }
-  }
 
   FixNH::init();
 }
@@ -93,9 +70,6 @@ void FixNHSphere::nve_v()
   double **torque = atom->torque;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
@@ -104,57 +78,18 @@ void FixNHSphere::nve_v()
 
   double dtfrotate = dtf / INERTIA;
   double dtirotate;
-  int itype;
 
   // update omega for all particles
   // d_omega/dt = torque / inertia
   // 4 cases depending on radius vs shape and rmass vs mass
 
-  if (radius) {
+  for (int i = 0; i < nlocal; i++)
-    if (rmass) {
+    if (mask[i] & groupbit) {
-      for (int i = 0; i < nlocal; i++) {    
+      dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
-	if (mask[i] & groupbit) {
+      omega[i][0] += dtirotate*torque[i][0];
-	  dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+      omega[i][1] += dtirotate*torque[i][1];
-	  omega[i][0] += dtirotate*torque[i][0];
+      omega[i][2] += dtirotate*torque[i][2];
-	  omega[i][1] += dtirotate*torque[i][1];
-	  omega[i][2] += dtirotate*torque[i][2];
-	}
-      }
-    } else {
-      for (int i = 0; i < nlocal; i++) {    
-	if (mask[i] & groupbit) {
-	  dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]);
-	  omega[i][0] += dtirotate*torque[i][0];
-	  omega[i][1] += dtirotate*torque[i][1];
-	  omega[i][2] += dtirotate*torque[i][2];
-	}
-      }
     }
-
-  } else {
-    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {    
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
-	  omega[i][0] += dtirotate*torque[i][0];
-	  omega[i][1] += dtirotate*torque[i][1];
-	  omega[i][2] += dtirotate*torque[i][2];
-	}
-      }
-    } else {
-      for (int i = 0; i < nlocal; i++) {    
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  dtirotate = dtfrotate / 
-	    (shape[itype][0]*shape[itype][0]*mass[itype]);
-	  omega[i][0] += dtirotate*torque[i][0];
-	  omega[i][1] += dtirotate*torque[i][1];
-	  omega[i][2] += dtirotate*torque[i][2];
-	}
-      }
-    }
-  }
 }
 
 /* ----------------------------------------------------------------------

src/fix_nve_sphere.cpp

@@ -53,10 +53,8 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) :
 
   // error checks
 
-  if (!atom->omega_flag || !atom->torque_flag)
+  if (!atom->sphere_flag)
-    error->all("Fix nve/sphere requires atom attributes omega, torque");
+    error->all("Fix nve/sphere requires atom style sphere");
-  if (!atom->radius_flag && !atom->avec->shape_type)
-    error->all("Fix nve/sphere requires atom attribute diameter or shape");
   if (extra == DIPOLE && !atom->mu_flag)
     error->all("Fix nve/sphere requires atom attribute mu");
 }
@@ -77,38 +75,17 @@ int FixNVESphere::setmask()
 
 void FixNVESphere::init()
 {
-  int i,itype;
+  // check that all particles are finite-size
-
-  // check that all particles are finite-size and spherical
   // no point particles allowed
 
-  if (atom->radius_flag) {
+  double *radius = atom->radius;
-    double *radius = atom->radius;
+  int *mask = atom->mask;
-    int *mask = atom->mask;
+  int nlocal = atom->nlocal;
-    int nlocal = atom->nlocal;
 
-    for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
+    if (mask[i] & groupbit)
-	if (radius[i] == 0.0)
+      if (radius[i] == 0.0)
-	  error->one("Fix nve/sphere requires extended particles");
+	error->one("Fix nve/sphere requires extended particles");
-      }
-
-  } else {
-    double **shape = atom->shape;
-    int *type = atom->type;
-    int *mask = atom->mask;
-    int nlocal = atom->nlocal;
-
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	itype = type[i];
-	if (shape[itype][0] == 0.0)
-	  error->one("Fix nve/sphere requires extended particles");
-	if (shape[itype][0] != shape[itype][1] || 
-	    shape[itype][0] != shape[itype][2])
-	  error->one("Fix nve/sphere requires spherical particle shapes");
-      }
-  }
 
   FixNVE::init();
 }
@@ -128,8 +105,6 @@ void FixNVESphere::initial_integrate(int vflag)
   double **torque = atom->torque;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  double **shape = atom->shape;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -141,86 +116,21 @@ void FixNVESphere::initial_integrate(int vflag)
 
   // update v,x,omega for all particles
   // d_omega/dt = torque / inertia
-  // 4 cases depending on radius vs shape and rmass vs mass
 
-  if (radius) {
+  for (int i = 0; i < nlocal; i++) {
-    if (rmass) {
+    if (mask[i] & groupbit) {
-      for (int i = 0; i < nlocal; i++) {
+      dtfm = dtf / rmass[i];
-	if (mask[i] & groupbit) {
+      v[i][0] += dtfm * f[i][0];
-	  dtfm = dtf / rmass[i];
+      v[i][1] += dtfm * f[i][1];
-	  v[i][0] += dtfm * f[i][0];
+      v[i][2] += dtfm * f[i][2];
-	  v[i][1] += dtfm * f[i][1];
+      x[i][0] += dtv * v[i][0];
-	  v[i][2] += dtfm * f[i][2];
+      x[i][1] += dtv * v[i][1];
-	  x[i][0] += dtv * v[i][0];
+      x[i][2] += dtv * v[i][2];
-	  x[i][1] += dtv * v[i][1];
+      
-	  x[i][2] += dtv * v[i][2];
+      dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
-	  
+      omega[i][0] += dtirotate * torque[i][0];
-	  dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+      omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][0] += dtirotate * torque[i][0];
+      omega[i][2] += dtirotate * torque[i][2];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
-
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  dtfm = dtf / mass[itype];
-	  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];
-	  
-	  dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]);
-	  omega[i][0] += dtirotate * torque[i][0];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
-    }
-
-  } else {
-    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  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];
-	  
-	  dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
-	  omega[i][0] += dtirotate * torque[i][0];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
-
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  dtfm = dtf / mass[itype];
-	  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];
-	  
-	  dtirotate = dtfrotate / 
-	    (shape[itype][0]*shape[itype][0]*mass[itype]);
-	  omega[i][0] += dtirotate * torque[i][0];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
     }
   }
 
@@ -230,21 +140,18 @@ void FixNVESphere::initial_integrate(int vflag)
 
   if (extra == DIPOLE) {
     double **mu = atom->mu;
-    double *dipole = atom->dipole;
+    for (int i = 0; i < nlocal; i++)
-    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit)
-      if (mask[i] & groupbit) {
+	if (mu[i][3] > 0.0) {
-	if (dipole[type[i]] > 0.0) {
 	  g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
 	  g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
 	  g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
 	  msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
-	  scale = dipole[type[i]]/sqrt(msq);
+	  scale = mu[i][3]/sqrt(msq);
 	  mu[i][0] = g[0]*scale;
 	  mu[i][1] = g[1]*scale;
 	  mu[i][2] = g[2]*scale;
 	}
-      }
-    }
   }
 }
 
@@ -252,18 +159,14 @@ void FixNVESphere::initial_integrate(int vflag)
 
 void FixNVESphere::final_integrate()
 {
-  int itype;
   double dtfm,dtirotate;
 
   double **v = atom->v;
   double **f = atom->f;
   double **omega = atom->omega;
   double **torque = atom->torque;
-  double *mass = atom->mass;
   double *rmass = atom->rmass;
   double *radius = atom->radius;
-  double **shape = atom->shape;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
@@ -274,74 +177,17 @@ void FixNVESphere::final_integrate()
 
   // update v,omega for all particles
   // d_omega/dt = torque / inertia
-  // 4 cases depending on radius vs shape and rmass vs mass
 
-  if (radius) {
+  for (int i = 0; i < nlocal; i++)
-    if (rmass) {
+    if (mask[i] & groupbit) {
-      for (int i = 0; i < nlocal; i++) {
+      dtfm = dtf / rmass[i];
-	if (mask[i] & groupbit) {
+      v[i][0] += dtfm * f[i][0];
-	  dtfm = dtf / rmass[i];
+      v[i][1] += dtfm * f[i][1];
-	  v[i][0] += dtfm * f[i][0];
+      v[i][2] += dtfm * f[i][2];
-	  v[i][1] += dtfm * f[i][1];
+      
-	  v[i][2] += dtfm * f[i][2];
+      dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
-	  
+      omega[i][0] += dtirotate * torque[i][0];
-	  dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+      omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][0] += dtirotate * torque[i][0];
+      omega[i][2] += dtirotate * torque[i][2];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
-
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  dtfm = dtf / mass[itype];
-	  v[i][0] += dtfm * f[i][0];
-	  v[i][1] += dtfm * f[i][1];
-	  v[i][2] += dtfm * f[i][2];
-	  
-	  dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]);
-	  omega[i][0] += dtirotate * torque[i][0];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
     }
-
-  } else {
-    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  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];
-	  
-	  dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
-	  omega[i][0] += dtirotate * torque[i][0];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
-
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-	if (mask[i] & groupbit) {
-	  itype = type[i];
-	  dtfm = dtf / mass[itype];
-	  v[i][0] += dtfm * f[i][0];
-	  v[i][1] += dtfm * f[i][1];
-	  v[i][2] += dtfm * f[i][2];
-	  
-	  dtirotate = dtfrotate / 
-	    (shape[itype][0]*shape[itype][0]*mass[itype]);
-	  omega[i][0] += dtirotate * torque[i][0];
-	  omega[i][1] += dtirotate * torque[i][1];
-	  omega[i][2] += dtirotate * torque[i][2];
-	}
-      }
-    }
-  }
 }

src/fix_rigid.cpp

@@ -352,8 +352,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
 
   // bitmasks for properties of extended particles
 
-  INERTIA_SPHERE_RADIUS = 1;
+  INERTIA_POINT = 1;
-  INERTIA_SPHERE_SHAPE = 2;
+  INERTIA_SPHERE = 2;
   INERTIA_ELLIPSOID = 4;
   ORIENT_DIPOLE = 8;
   ORIENT_QUAT = 16;
@@ -462,30 +462,32 @@ void FixRigid::init()
     step_respa = ((Respa *) update->integrate)->step;
 
   // extended = 1 if any particle in a rigid body is finite size
+  //              or has a dipole moment
 
   extended = dorientflag = qorientflag = 0;
 
+  double **shape = atom->shape;
+  double **mu = atom->mu;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
-  double **shape = atom->shape;
-  double *dipole = atom->dipole;
   int *type = atom->type;
   int nlocal = atom->nlocal;
 
-  if (atom->radius_flag || atom->avec->shape_type) {
+  if (atom->radius_flag || atom->shape_flag) {
     int flag = 0;
     for (i = 0; i < nlocal; i++) {
       if (body[i] < 0) continue;
       if (radius && radius[i] > 0.0) flag = 1;
-      else if (shape && shape[type[i]][0] > 0.0) flag = 1;
+      if (shape && shape[i][0] > 0.0) flag = 1;
+      if (mu && mu[i][3] > 0.0) flag = 1;
     }
 
     MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
   }
 
   // grow extended arrays and set extended flags for each particle
-  // vorientflag = 1 if any particles store dipole orientation
+  // dorientflag = 1 if any particles store dipole orientation
   // qorientflag = 1 if any particles store quat orientation
 
   if (extended) {
@@ -497,39 +499,23 @@ void FixRigid::init()
       eflags[i] = 0;
       if (body[i] < 0) continue;
 
-      // set INERTIA if radius or shape > 0.0
+      // set INERTIA to POINT or SPHERE or ELLIPSOID
 
-      if (radius) {
+      if (radius && radius[i] > 0.0) {
-	if (radius[i] > 0.0) eflags[i] |= INERTIA_SPHERE_RADIUS;
+	eflags[i] |= INERTIA_SPHERE;
-      } else if (shape) {
+	eflags[i] |= OMEGA;
-	if (shape[type[i]][0] > 0.0) {
+	eflags[i] |= TORQUE;
-	  if (shape[type[i]][0] == shape[type[i]][1] &&
+      } else if (shape && shape[i][0] > 0.0) {
-	      shape[type[i]][0] == shape[type[i]][2])
+	eflags[i] |= INERTIA_ELLIPSOID;
-	    eflags[i] |= INERTIA_SPHERE_SHAPE;
+	eflags[i] |= ORIENT_QUAT;
-	  else eflags[i] |= INERTIA_ELLIPSOID;
+	eflags[i] |= ANGMOM;
-	}
+	eflags[i] |= TORQUE;
-      }
+      } else eflags[i] |= INERTIA_POINT;
 
-      // other flags only set if particle is finite size
+      // set DIPOLE if atom->mu and mu[3] > 0.0
-      // set DIPOLE if atom->mu and atom->dipole exist and dipole[itype] > 0.0
-      // set QUAT if atom->quat exists (could be ellipsoid or sphere)
-      // set TORQUE if atom->torque exists
-      // set exactly one of OMEGA and ANGMOM so particle contributes once
-      // set OMEGA if either radius or rmass exists
-      // set ANGMOM if shape and mass exist
-      // set OMEGA if atom->angmom doesn't exist
 
-      if (eflags[i] == 0) continue;
+      if (atom->mu_flag && mu[i][3] > 0.0)
-
-      if (atom->mu_flag && dipole && dipole[type[i]] > 0.0)
 	eflags[i] |= ORIENT_DIPOLE;
-      if (atom->quat_flag) eflags[i] |= ORIENT_QUAT;
-      if (atom->torque_flag) eflags[i] |= TORQUE;
-      if ((radius || rmass) && atom->omega_flag) eflags[i] |= OMEGA;
-      else if (shape && mass && atom->angmom_flag) eflags[i] |= ANGMOM;
-      else if (atom->omega_flag) eflags[i] |= OMEGA;
-      else error->one("Could not set finite-size particle attribute "
-		      "in fix rigid");
     }
   }
 
@@ -638,7 +624,7 @@ void FixRigid::init()
     sum[ibody][5] -= massone * dx*dz;
   }
 
-  // extended particles contribute extra terms to moments of inertia
+  // extended particles may contribute extra terms to moments of inertia
 
   if (extended) {
     double ex[3],ey[3],ez[3],idiag[3];
@@ -655,26 +641,20 @@ void FixRigid::init()
       if (rmass) massone = rmass[i];
       else massone = mass[itype];
 
-      if (eflags[i] & INERTIA_SPHERE_RADIUS) {
+      if (eflags[i] & INERTIA_SPHERE) {
 	sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
       }
-      if (eflags[i] & INERTIA_SPHERE_SHAPE) {
-	rad = shape[type[i]][0];
-	sum[ibody][0] += 0.4 * massone * rad*rad;
-	sum[ibody][1] += 0.4 * massone * rad*rad;
-	sum[ibody][2] += 0.4 * massone * rad*rad;
-      }
       if (eflags[i] & INERTIA_ELLIPSOID) {
 	MathExtra::quat_to_mat(atom->quat[i],p);
 	MathExtra::quat_to_mat_trans(atom->quat[i],ptrans);
 	idiag[0] = 0.2*massone *
-	  (shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]);
+	  (shape[i][1]*shape[i][1]+shape[itype][2]*shape[i][2]);
 	idiag[1] = 0.2*massone *
-	  (shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]);
+	  (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]);
 	idiag[2] = 0.2*massone *
-	  (shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][1]);
+	  (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]);
 	MathExtra::diag_times3(idiag,ptrans,itemp);
 	MathExtra::times3(p,itemp,ispace);
 	sum[ibody][0] += ispace[0][0];
@@ -800,7 +780,6 @@ void FixRigid::init()
     
     if (extended) {
       double **mu = atom->mu;
-      double *dipole = atom->dipole;
       int *type = atom->type;
 
       if (eflags[i] & ORIENT_DIPOLE) {
@@ -810,7 +789,7 @@ void FixRigid::init()
 	  mu[i][1]*ey_space[ibody][1] + mu[i][2]*ey_space[ibody][2];
 	dorient[i][2] = mu[i][0]*ez_space[ibody][0] + 
 	  mu[i][1]*ez_space[ibody][1] + mu[i][2]*ez_space[ibody][2];
-	MathExtra::snormalize3(dipole[type[i]],dorient[i],dorient[i]);
+	MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]);
       } else if (dorientflag)
 	dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
 
@@ -827,7 +806,7 @@ void FixRigid::init()
   // recompute moments of inertia around new axes
   // 3 diagonal moments should equal principal moments
   // 3 off-diagonal moments should be 0.0
-  // extended particles contribute extra terms to moments of inertia
+  // extended particles may contribute extra terms to moments of inertia
   
   for (ibody = 0; ibody < nbody; ibody++)
     for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
@@ -864,26 +843,20 @@ void FixRigid::init()
       if (rmass) massone = rmass[i];
       else massone = mass[itype];
 
-      if (eflags[i] & INERTIA_SPHERE_RADIUS) {
+      if (eflags[i] & INERTIA_SPHERE) {
 	sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
       }
-      if (eflags[i] & INERTIA_SPHERE_SHAPE) {
-	rad = shape[type[i]][0];
-	sum[ibody][0] += 0.4 * massone * rad*rad;
-	sum[ibody][1] += 0.4 * massone * rad*rad;
-	sum[ibody][2] += 0.4 * massone * rad*rad;
-      }
       if (eflags[i] & INERTIA_ELLIPSOID) {
 	MathExtra::quat_to_mat(qorient[i],p);
 	MathExtra::quat_to_mat_trans(qorient[i],ptrans);
 	idiag[0] = 0.2*massone *
-	  (shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]);
+	  (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]);
 	idiag[1] = 0.2*massone *
-	  (shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]);
+	  (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]);
 	idiag[2] = 0.2*massone *
-	  (shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][1]);
+	  (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]);
 	MathExtra::diag_times3(idiag,ptrans,itemp);
 	MathExtra::times3(p,itemp,ispace);
 	sum[ibody][0] += ispace[0][0];
@@ -1032,7 +1005,6 @@ void FixRigid::setup(int vflag)
     double **angmom_one = atom->angmom;
     double **torque_one = atom->torque;
     double *radius = atom->radius;
-    double **shape = atom->shape;
 
     for (i = 0; i < nlocal; i++) {
       if (body[i] < 0) continue;
@@ -1041,8 +1013,7 @@ void FixRigid::setup(int vflag)
       if (eflags[i] & OMEGA) {
 	if (rmass) massone = rmass[i];
 	else massone = mass[type[i]];
-	if (radius) radone = radius[i];
+	radone = radius[i];
-	else radone = shape[type[i]][0];
 	sum[ibody][0] += 0.4 * massone * radone*radone * omega_one[i][0];
 	sum[ibody][1] += 0.4 * massone * radone*radone * omega_one[i][1];
 	sum[ibody][2] += 0.4 * massone * radone*radone * omega_one[i][2];
@@ -1965,8 +1936,8 @@ void FixRigid::set_xv()
   if (extended) {
     double **omega_one = atom->omega;
     double **angmom_one = atom->angmom;
-    double *dipole = atom->dipole;
     double **shape = atom->shape;
+    double **mu = atom->mu;
 
     for (int i = 0; i < nlocal; i++) {
       if (body[i] < 0) continue;
@@ -1974,8 +1945,8 @@ void FixRigid::set_xv()
       
       if (eflags[i] & ORIENT_DIPOLE) {
 	MathExtra::quat_to_mat(quat[ibody],p);
-	MathExtra::times_column3(p,dorient[i],atom->mu[i]);
+	MathExtra::times_column3(p,dorient[i],mu[i]);
-	MathExtra::snormalize3(dipole[type[i]],atom->mu[i],atom->mu[i]);
+	MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
       }
       if (eflags[i] & ORIENT_QUAT) {
 	quatquat(quat[ibody],qorient[i],atom->quat[i]);
@@ -1987,13 +1958,12 @@ void FixRigid::set_xv()
 	omega_one[i][2] = omega[ibody][2];
       }
       if (eflags[i] & ANGMOM) {
-	itype = type[i];
+	ione[0] = 0.2*rmass[i] *
-	ione[0] = 0.2*mass[itype] *
+	  (shape[i][1]*shape[i][1] + shape[i][2]*shape[i][2]);
-	  (shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]);
+	ione[1] = 0.2*rmass[i] *
-	ione[1] = 0.2*mass[itype] *
+	  (shape[i][0]*shape[i][0] + shape[i][2]*shape[i][2]);
-	  (shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]);
+	ione[2] = 0.2*rmass[i] * 
-	ione[2] = 0.2*mass[itype] * 
+	  (shape[i][0]*shape[i][0] + shape[i][1]*shape[i][1]);
-	  (shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]);
 	exyz_from_q(atom->quat[i],exone,eyone,ezone);
 	angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
       }
@@ -2117,13 +2087,12 @@ void FixRigid::set_v()
 	omega_one[i][2] = omega[ibody][2];
       }
       if (eflags[i] & ANGMOM) {
-	itype = type[i];
+	ione[0] = 0.2*rmass[i] *
-	ione[0] = 0.2*mass[itype] *
+	  (shape[i][1]*shape[i][1] + shape[i][2]*shape[i][2]);
-	  (shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]);
+	ione[1] = 0.2*rmass[i] *
-	ione[1] = 0.2*mass[itype] *
+	  (shape[i][0]*shape[i][0] + shape[i][2]*shape[i][2]);
-	  (shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]);
+	ione[2] = 0.2*rmass[i] * 
-	ione[2] = 0.2*mass[itype] * 
+	  (shape[i][0]*shape[i][0] + shape[i][1]*shape[i][1]);
-	  (shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]);
 	exyz_from_q(atom->quat[i],exone,eyone,ezone);
 	angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
       }

src/fix_rigid.h

@@ -96,7 +96,7 @@ class FixRigid : public Fix {
   int p_chain;
 
                             // bitmasks for eflags
-  int INERTIA_SPHERE_RADIUS,INERTIA_SPHERE_SHAPE,INERTIA_ELLIPSOID;
+  int INERTIA_POINT,INERTIA_SPHERE,INERTIA_ELLIPSOID;
   int ORIENT_DIPOLE,ORIENT_QUAT;
   int OMEGA,ANGMOM,TORQUE;
 

src/fix_wall_region.cpp

@@ -98,32 +98,20 @@ void FixWallRegion::init()
     error->all("Region ID for fix wall/region does not exist");
 
   // error checks for style COLLOID
-  // insure all particle shapes are spherical
-  // can be polydisperse
   // insure all particles in group are extended particles
 
   if (style == COLLOID) {
-    if (!atom->avec->shape_type)
+    if (!atom->sphere_flag) 
-      error->all("Fix wall/region colloid requires atom attribute shape");
+      error->all("Fix wall/region colloid requires atom style sphere");
-    if (atom->radius_flag)
-      error->all("Fix wall/region colloid cannot be used with "
-		 "atom attribute diameter");
 
-    for (int i = 1; i <= atom->ntypes; i++)
+    double *radius = atom->radius;
-      if ((atom->shape[i][0] != atom->shape[i][1]) || 
-	  (atom->shape[i][0] != atom->shape[i][2]) ||
-	  (atom->shape[i][1] != atom->shape[i][2]))
-	error->all("Fix wall/region colloid requires spherical particles");
-
-    double **shape = atom->shape;
-    int *type = atom->type;
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
     
     int flag = 0;
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit)
-	if (atom->shape[type[i]][0] == 0.0) flag = 1;
+	if (radius[i] == 0.0) flag = 1;
     
     int flagall;
     MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
@@ -197,7 +185,7 @@ void FixWallRegion::post_force(int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double **shape = atom->shape;
+  double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -207,7 +195,7 @@ void FixWallRegion::post_force(int vflag)
 
   // region->match() insures particle is in region or on surface, else error
   // if returned contact dist r = 0, is on surface, also an error
-  // in COLLOID case, r <= shape radius is an error
+  // in COLLOID case, r <= radius is an error
   
   for (i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
@@ -215,7 +203,7 @@ void FixWallRegion::post_force(int vflag)
 	onflag = 1;
 	continue;
       }
-      if (style == COLLOID) tooclose = shape[type[i]][0];
+      if (style == COLLOID) tooclose = radius[i];
       else tooclose = 0.0;
 
       n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);
@@ -228,8 +216,7 @@ void FixWallRegion::post_force(int vflag)
 
 	if (style == LJ93) lj93(region->contact[m].r);
 	else if (style == LJ126) lj126(region->contact[m].r);
-	else if (style == COLLOID) 
+	else if (style == COLLOID) colloid(region->contact[m].r,radius[i]);
-	  colloid(region->contact[m].r,shape[type[i]][0]);
 	else harmonic(region->contact[m].r);
 
 	ewall[0] += eng;

src/input.cpp

@@ -431,7 +431,6 @@ int Input::execute_command()
   else if (!strcmp(command,"dihedral_coeff")) dihedral_coeff();
   else if (!strcmp(command,"dihedral_style")) dihedral_style();
   else if (!strcmp(command,"dimension")) dimension();
-  else if (!strcmp(command,"dipole")) dipole();
   else if (!strcmp(command,"dump")) dump();
   else if (!strcmp(command,"dump_modify")) dump_modify();
   else if (!strcmp(command,"fix")) fix();
@@ -457,7 +456,6 @@ int Input::execute_command()
   else if (!strcmp(command,"reset_timestep")) reset_timestep();
   else if (!strcmp(command,"restart")) restart();
   else if (!strcmp(command,"run_style")) run_style();
-  else if (!strcmp(command,"shape")) shape();
   else if (!strcmp(command,"special_bonds")) special_bonds();
   else if (!strcmp(command,"thermo")) thermo();
   else if (!strcmp(command,"thermo_modify")) thermo_modify();
@@ -948,16 +946,6 @@ void Input::dimension()
 
 /* ---------------------------------------------------------------------- */
 
-void Input::dipole()
-{
-  if (narg != 2) error->all("Illegal dipole command");
-  if (domain->box_exist == 0)
-    error->all("Dipole command before simulation box is defined");
-  atom->set_dipole(narg,arg);
-}
-
-/* ---------------------------------------------------------------------- */
-
 void Input::dump()
 {
   output->add_dump(narg,arg);
@@ -1208,16 +1196,6 @@ void Input::run_style()
 
 /* ---------------------------------------------------------------------- */
 
-void Input::shape()
-{
-  if (narg != 4) error->all("Illegal shape command");
-  if (domain->box_exist == 0)
-    error->all("Shape command before simulation box is defined");
-  atom->set_shape(narg,arg);
-}
-
-/* ---------------------------------------------------------------------- */
-
 void Input::special_bonds()
 {
   // store 1-3,1-4 and dihedral/extra flag values before change

src/input.h

@@ -75,7 +75,6 @@ class Input : protected Pointers {
   void dihedral_coeff();
   void dihedral_style();
   void dimension();
-  void dipole();
   void dump();
   void dump_modify();
   void fix();
@@ -101,7 +100,6 @@ class Input : protected Pointers {
   void reset_timestep();
   void restart();
   void run_style();
-  void shape();
   void special_bonds();
   void thermo();
   void thermo_modify();

src/lammps.cpp

@@ -338,7 +338,7 @@ void LAMMPS::init()
                          //   atom_vec init uses deform_vremap
   modify->init();        // modify must come after update, force, atom, domain
   neighbor->init();      // neighbor must come after force, modify
-  comm->init();          // comm must come after force, modify, neighbor
+  comm->init();          // comm must come after force, modify, neighbor, atom
   output->init();        // output must come after domain, force, modify
   timer->init();
 }

src/pair_yukawa.cpp

@@ -38,6 +38,7 @@ PairYukawa::~PairYukawa()
     memory->destroy(setflag);
     memory->destroy(cutsq);
 
+    memory->destroy(rad);
     memory->destroy(cut);
     memory->destroy(a);
     memory->destroy(offset);
@@ -139,6 +140,7 @@ void PairYukawa::allocate()
 
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
 
+  memory->create(rad,n+1,"pair:rad");
   memory->create(cut,n+1,n+1,"pair:cut");
   memory->create(a,n+1,n+1,"pair:a");
   memory->create(offset,n+1,n+1,"pair:offset");

src/pair_yukawa.h

@@ -41,6 +41,7 @@ class PairYukawa : public Pair {
  protected:
   double cut_global;
   double kappa;
+  double *rad;
   double **cut,**a,**offset;
 
   void allocate();

src/read_data.cpp

@@ -163,10 +163,6 @@ void ReadData::command(int narg, char **arg)
 
     } else if (strcmp(keyword,"Masses") == 0) {
       mass();
-    } else if (strcmp(keyword,"Shapes") == 0) {
-      shape();
-    } else if (strcmp(keyword,"Dipoles") == 0) {
-      dipole();
     } else if (strcmp(keyword,"Pair Coeffs") == 0) {
       if (force->pair == NULL) 
 	error->all("Must define pair_style before Pair Coeffs");
@@ -292,8 +288,7 @@ void ReadData::header(int flag)
   char *ptr;
 
   char *section_keywords[NSECTIONS] = 
-    {"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers",
+    {"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers","Masses",
-     "Masses","Shapes","Dipoles",
      "Pair Coeffs","Bond Coeffs","Angle Coeffs",
      "Dihedral Coeffs","Improper Coeffs",
      "BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs",
@@ -761,64 +756,6 @@ void ReadData::mass()
 
 /* ---------------------------------------------------------------------- */
 
-void ReadData::shape()
-{
-  int i,m;
-  char *buf = new char[atom->ntypes*MAXLINE];
-  char *original = buf;
-
-  if (me == 0) {
-    char *eof;
-    m = 0;
-    for (i = 0; i < atom->ntypes; i++) {
-      eof = fgets(&buf[m],MAXLINE,fp);
-      if (eof == NULL) error->one("Unexpected end of data file");
-      m += strlen(&buf[m]);
-      buf[m-1] = '\0';
-    }
-  }
-
-  MPI_Bcast(&m,1,MPI_INT,0,world);
-  MPI_Bcast(buf,m,MPI_CHAR,0,world);
-
-  for (i = 0; i < atom->ntypes; i++) {
-    atom->set_shape(buf);
-    buf += strlen(buf) + 1;
-  }
-  delete [] original;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ReadData::dipole()
-{
-  int i,m;
-  char *buf = new char[atom->ntypes*MAXLINE];
-  char *original = buf;
-
-  if (me == 0) {
-    char *eof;
-    m = 0;
-    for (i = 0; i < atom->ntypes; i++) {
-      eof = fgets(&buf[m],MAXLINE,fp);
-      if (eof == NULL) error->one("Unexpected end of data file");
-      m += strlen(&buf[m]);
-      buf[m-1] = '\0';
-    }
-  }
-
-  MPI_Bcast(&m,1,MPI_INT,0,world);
-  MPI_Bcast(buf,m,MPI_CHAR,0,world);
-
-  for (i = 0; i < atom->ntypes; i++) {
-    atom->set_dipole(buf);
-    buf += strlen(buf) + 1;
-  }
-  delete [] original;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void ReadData::paircoeffs()
 {
   int i,m;
@@ -1007,7 +944,6 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
   while (strlen(keyword)) {
 
     if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes);
-    else if (strcmp(keyword,"Dipoles") == 0) skip_lines(atom->ntypes);
     else if (strcmp(keyword,"Atoms") == 0) skip_lines(natoms);
     else if (strcmp(keyword,"Velocities") == 0) skip_lines(natoms);
 

src/read_data.h

@@ -54,7 +54,6 @@ class ReadData : protected Pointers {
   void impropers();
 
   void mass();
-  void shape();
   void dipole();
 
   void paircoeffs();

src/read_restart.cpp

@@ -56,7 +56,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT,
        SPECIAL_LJ_1,SPECIAL_LJ_2,SPECIAL_LJ_3,
        SPECIAL_COUL_1,SPECIAL_COUL_2,SPECIAL_COUL_3,
        XY,XZ,YZ};
-enum{MASS,SHAPE,DIPOLE};
+enum{MASS};
 enum{PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
 
 #define LB_FACTOR 1.1
@@ -700,21 +700,6 @@ void ReadRestart::type_arrays()
       atom->set_mass(mass);
       delete [] mass;
 
-    } else if (flag == SHAPE) {
-      double **shape;
-      memory->create(shape,atom->ntypes+1,3,"restart:shape");
-      if (me == 0) fread(&shape[1][0],sizeof(double),atom->ntypes*3,fp);
-      MPI_Bcast(&shape[1][0],atom->ntypes*3,MPI_DOUBLE,0,world);
-      atom->set_shape(shape);
-      memory->destroy(shape);
-
-    } else if (flag == DIPOLE) {
-      double *dipole = new double[atom->ntypes+1];
-      if (me == 0) fread(&dipole[1],sizeof(double),atom->ntypes,fp);
-      MPI_Bcast(&dipole[1],atom->ntypes,MPI_DOUBLE,0,world);
-      atom->set_dipole(dipole);
-      delete [] dipole;
-
     } else error->all("Invalid flag in type arrays section of restart file");
 
     flag = read_int();

src/replicate.cpp

@@ -208,25 +208,6 @@ void Replicate::command(int narg, char **arg)
     }
   }
 
-  if (atom->shape) {
-    for (int itype = 1; itype <= atom->ntypes; itype++) {
-      atom->shape_setflag[itype] = old->shape_setflag[itype];
-      if (atom->shape_setflag[itype]) {
-	atom->shape[itype][0] = old->shape[itype][0];
-	atom->shape[itype][1] = old->shape[itype][1];
-	atom->shape[itype][2] = old->shape[itype][2];
-      }
-    }
-  }
-
-  if (atom->dipole) {
-    for (int itype = 1; itype <= atom->ntypes; itype++) {
-      atom->dipole_setflag[itype] = old->dipole_setflag[itype];
-      if (atom->dipole_setflag[itype]) 
-	atom->dipole[itype] = old->dipole[itype];
-    }
-  }
-
   // set bounds for my proc
   // if periodic and I am lo/hi proc, adjust bounds by EPSILON
   // insures all replicated atoms will be owned even with round-off

src/set.cpp

@@ -31,12 +31,12 @@
 
 using namespace LAMMPS_NS;
 
-enum{ATOM,GROUP,REGION};
+enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
-enum{TYPE,TYPE_FRACTION,MOLECULE,
+enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,
-       X,Y,Z,CHARGE,
+     DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,
-       DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,
+     DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER};
-       DIAMETER,DENSITY,VOLUME,IMAGE,
+
-       BOND,ANGLE,DIHEDRAL,IMPROPER};
+#define BIG INT_MAX
 
 /* ---------------------------------------------------------------------- */
 
@@ -55,11 +55,13 @@ void Set::command(int narg, char **arg)
     error->all("Set command with no atoms existing");
   if (narg < 3) error->all("Illegal set command");
 
-  // style and ID
+  // style and ID info
 
-  if (strcmp(arg[0],"atom") == 0) style = ATOM;
+  if (strcmp(arg[0],"atom") == 0) style = ATOM_SELECT;
-  else if (strcmp(arg[0],"group") == 0) style = GROUP;
+  else if (strcmp(arg[0],"mol") == 0) style = MOL_SELECT;
-  else if (strcmp(arg[0],"region") == 0) style = REGION;
+  else if (strcmp(arg[0],"type") == 0) style = TYPE_SELECT;
+  else if (strcmp(arg[0],"group") == 0) style = GROUP_SELECT;
+  else if (strcmp(arg[0],"region") == 0) style = REGION_SELECT;
   else error->all("Illegal set command");
 
   int n = strlen(arg[1]) + 1;
@@ -91,6 +93,8 @@ void Set::command(int narg, char **arg)
       newtype = atoi(arg[iarg+1]);
       fraction = atof(arg[iarg+2]);
       ivalue = atoi(arg[iarg+3]);
+      if (newtype <= 0 || newtype > atom->ntypes)
+	error->all("Invalid value in set command");
       if (fraction < 0.0 || fraction > 1.0) 
 	error->all("Invalid value in set command");
       if (ivalue <= 0) error->all("Invalid random number seed in set command");
@@ -125,23 +129,48 @@ void Set::command(int narg, char **arg)
 	error->all("Cannot set this attribute for this atom style");
       set(CHARGE);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"mass") == 0) {
+      if (iarg+2 > narg) error->all("Illegal set command");
+      dvalue = atof(arg[iarg+1]);
+      if (!atom->rmass_flag)
+	error->all("Cannot set this attribute for this atom style");
+      if (dvalue <= 0.0) error->all("Invalid mass in set command");
+      set(MASS);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"shape") == 0) {
+      if (iarg+4 > narg) error->all("Illegal set command");
+      xvalue = atof(arg[iarg+1]);
+      yvalue = atof(arg[iarg+2]);
+      zvalue = atof(arg[iarg+3]);
+      if (!atom->shape_flag)
+	error->all("Cannot set this attribute for this atom style");
+      if (xvalue < 0.0 || yvalue < 0.0 || zvalue < 0.0)
+	error->all("Invalid shape in set command");
+      if (xvalue > 0.0 || yvalue > 0.0 || zvalue > 0.0) {
+	if (xvalue == 0.0 || yvalue == 0.0 || zvalue == 0.0)
+	  error->one("Invalid shape in set command");
+      }
+      set(SHAPE);
+      iarg += 4;
     } else if (strcmp(arg[iarg],"dipole") == 0) {
       if (iarg+4 > narg) error->all("Illegal set command");
       xvalue = atof(arg[iarg+1]);
       yvalue = atof(arg[iarg+2]);
       zvalue = atof(arg[iarg+3]);
-      if (!atom->mu_flag || atom->dipole == NULL)
+      if (!atom->mu_flag)
 	error->all("Cannot set this attribute for this atom style");
       set(DIPOLE);
       iarg += 4;
     } else if (strcmp(arg[iarg],"dipole/random") == 0) {
-      if (iarg+2 > narg) error->all("Illegal set command");
+      if (iarg+3 > narg) error->all("Illegal set command");
       ivalue = atoi(arg[iarg+1]);
-      if (!atom->mu_flag || atom->shape == NULL)
+      dvalue = atof(arg[iarg+2]);
+      if (!atom->mu_flag)
 	error->all("Cannot set this attribute for this atom style");
       if (ivalue <= 0) error->all("Invalid random number seed in set command");
+      if (dvalue <= 0.0) error->all("Invalid dipole length in set command");
       setrandom(DIPOLE_RANDOM);
-      iarg += 2;
+      iarg += 3;
     } else if (strcmp(arg[iarg],"quat") == 0) {
       if (iarg+5 > narg) error->all("Illegal set command");
       xvalue = atof(arg[iarg+1]);
@@ -165,12 +194,13 @@ void Set::command(int narg, char **arg)
       dvalue = atof(arg[iarg+1]);
       if (!atom->radius_flag)
         error->all("Cannot set this attribute for this atom style");
+      if (dvalue < 0.0) error->all("Invalid diameter in set command");
       set(DIAMETER);
       iarg += 2;
     } else if (strcmp(arg[iarg],"density") == 0) {
       if (iarg+2 > narg) error->all("Illegal set command");
       dvalue = atof(arg[iarg+1]);
-      if (!atom->density_flag)
+      if (!atom->rmass_flag)
         error->all("Cannot set this attribute for this atom style");
       set(DENSITY);
       iarg += 2;
@@ -261,7 +291,7 @@ void Set::command(int narg, char **arg)
 }
 
 /* ----------------------------------------------------------------------
-   select atoms according to ATOM, GROUP, REGION style
+   select atoms according to ATOM, MOLECULE, TYPE, GROUP, REGION style
    n = nlocal or nlocal+nghost depending on keyword
 ------------------------------------------------------------------------- */
 
@@ -269,18 +299,38 @@ void Set::selection(int n)
 {
   delete [] select;
   select = new int[n];
+  int nlo,nhi;
 
-  if (style == ATOM) {
+  if (style == ATOM_SELECT) {
     if (atom->tag_enable == 0)
       error->all("Cannot use set atom with no atom IDs defined");
-    int idatom = atoi(id);
+    force->bounds(id,BIG,nlo,nhi);
 
     int *tag = atom->tag;
     for (int i = 0; i < n; i++)
-      if (idatom == tag[i]) select[i] = 1;
+      if (tag[i] >= nlo && tag[i] <= nhi) select[i] = 1;
       else select[i] = 0;
 
-  } else if (style == GROUP) {
+  } else if (style == MOL_SELECT) {
+    if (atom->molecule_flag == 0)
+      error->all("Cannot use set mol with no molecule IDs defined");
+    if (strcmp(id,"0") == 0) nlo = nhi = 0;
+    else force->bounds(id,BIG,nlo,nhi);
+
+    int *molecule = atom->molecule;
+    for (int i = 0; i < n; i++)
+      if (molecule[i] >= nlo && molecule[i] <= nhi) select[i] = 1;
+      else select[i] = 0;
+
+  } else if (style == TYPE_SELECT) {
+    force->bounds(id,atom->ntypes,nlo,nhi);
+
+    int *type = atom->type;
+    for (int i = 0; i < n; i++)
+      if (type[i] >= nlo && type[i] <= nhi) select[i] = 1;
+      else select[i] = 0;
+
+  } else if (style == GROUP_SELECT) {
     int igroup = group->find(id);
     if (igroup == -1) error->all("Could not find set group ID");
     int groupbit = group->bitmask[igroup];
@@ -290,7 +340,7 @@ void Set::selection(int n)
       if (mask[i] & groupbit) select[i] = 1;
       else select[i] = 0;
 
-  } else {
+  } else if (style == REGION_SELECT) {
     int iregion = domain->find_region(id);
     if (iregion == -1) error->all("Set region ID does not exist");
 
@@ -308,83 +358,81 @@ void Set::selection(int n)
 
 void Set::set(int keyword)
 {
-  if (keyword == DIPOLE) atom->check_dipole();
-
   selection(atom->nlocal);
 
   int nlocal = atom->nlocal;
-  for (int i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++) {
-    if (select[i]) {
+    if (!select[i]) continue;
-      if (keyword == TYPE) atom->type[i] = ivalue;
-      else if (keyword == MOLECULE) atom->molecule[i] = ivalue;
-      else if (keyword == X) atom->x[i][0] = dvalue;
-      else if (keyword == Y) atom->x[i][1] = dvalue;
-      else if (keyword == Z) atom->x[i][2] = dvalue;
-      else if (keyword == CHARGE) atom->q[i] = dvalue;
 
-      // set radius from diameter
+    if (keyword == TYPE) atom->type[i] = ivalue;
-      // set rmass if both rmass and density are defined
+    else if (keyword == MOLECULE) atom->molecule[i] = ivalue;
+    else if (keyword == X) atom->x[i][0] = dvalue;
+    else if (keyword == Y) atom->x[i][1] = dvalue;
+    else if (keyword == Z) atom->x[i][2] = dvalue;
+    else if (keyword == CHARGE) atom->q[i] = dvalue;
+    else if (keyword == MASS) atom->rmass[i] = dvalue;
+    else if (keyword == DIAMETER) atom->radius[i] = 0.5 * dvalue;
+    else if (keyword == VOLUME) atom->vfrac[i] = dvalue;
+    
+    // set shape
 
-      else if (keyword == DIAMETER) {
+    else if (keyword == SHAPE) {
-	atom->radius[i] = 0.5 * dvalue;
+      double **shape = atom->shape;
-	if (atom->rmass_flag && atom->density_flag) 
+      shape[i][0] = 0.5 * xvalue;
-	  atom->rmass[i] = 4.0*PI/3.0 * 
+      shape[i][1] = 0.5 * yvalue;
-	    atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i];
+      shape[i][2] = 0.5 * zvalue;
+      
+    // set rmass via density
+    // if radius > 0.0, assume sphere
+    // if shape > 0.0, assume ellipsoid
+    // else set rmass to density directly
 
-      // set density
+    } else if (keyword == DENSITY) {
-      // set rmass (granular) if both rmass and radius are defined
+      if (atom->radius_flag && atom->radius[i] > 0.0)
-      // set rmass (peri) if both rmass and vfrac are defined
+	atom->rmass[i] = 4.0*PI/3.0 * 
-	
+	  atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
-      } else if (keyword == DENSITY) {
+      else if (atom->shape_flag && atom->shape[i][0] > 0.0)
-	atom->density[i] = dvalue;
+	atom->rmass[i] = 4.0*PI/3.0 * 
-	if (atom->rmass_flag && atom->radius_flag) 
+	  atom->shape[i][0]*atom->shape[i][1]*atom->shape[i][2] * dvalue;
-	  atom->rmass[i] = 4.0*PI/3.0 * 
+      else atom->rmass[i] = dvalue;
-	    atom->radius[i]*atom->radius[i]*atom->radius[i] * 
-	    atom->density[i];
-	else if (atom->rmass_flag && atom->vfrac_flag)
-	  atom->rmass[i] = dvalue;
 
-      } else if (keyword == VOLUME) atom->vfrac[i] = dvalue;
+    // reset any or all of 3 image flags
-
+      
-      // reset any or all of 3 image flags
+    } else if (keyword == IMAGE) {
-
+      int xbox = (atom->image[i] & 1023) - 512;
-      else if (keyword == IMAGE) {
+      int ybox = (atom->image[i] >> 10 & 1023) - 512;
-	int xbox = (atom->image[i] & 1023) - 512;
+      int zbox = (atom->image[i] >> 20) - 512;
-	int ybox = (atom->image[i] >> 10 & 1023) - 512;
+      if (ximageflag) xbox = ximage;
-	int zbox = (atom->image[i] >> 20) - 512;
+      if (yimageflag) ybox = yimage;
-	if (ximageflag) xbox = ximage;
+      if (zimageflag) zbox = zimage;
-	if (yimageflag) ybox = yimage;
+      atom->image[i] = ((zbox + 512 & 1023) << 20) |
-	if (zimageflag) zbox = zimage;
+	((ybox + 512 & 1023) << 10) |	(xbox + 512 & 1023);
-	atom->image[i] = ((zbox + 512 & 1023) << 20) |
+      
-	  ((ybox + 512 & 1023) << 10) |	(xbox + 512 & 1023);
+    // set dipole moment
-
+      
-      } else if (keyword == DIPOLE) {
+    } else if (keyword == DIPOLE) {
-	if (atom->dipole[atom->type[i]] > 0.0) {
+      double **mu = atom->mu;
-	  double **mu = atom->mu;
+      mu[i][0] = xvalue;
-	  mu[i][0] = xvalue;
+      mu[i][1] = yvalue;
-	  mu[i][1] = yvalue;
+      mu[i][2] = zvalue;
-	  mu[i][2] = zvalue;
+      mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + 
-	  double lensq = 
+		      mu[i][2]*mu[i][2]);
-	    mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
+      
-	  double scale = atom->dipole[atom->type[i]]/sqrt(lensq);
+    // set quaternion orientation
-	  mu[i][0] *= scale;
+      
-	  mu[i][1] *= scale;
+    } else if (keyword == QUAT) {
-	  mu[i][2] *= scale;
+      double PI = 4.0*atan(1.0);
-	}
+      double theta2 = 0.5 * PI * wvalue/180.0;
-
+      double sintheta2 = sin(theta2);
-      } else if (keyword == QUAT) {
+      double **quat = atom->quat;
-	double PI = 4.0*atan(1.0);
+      quat[i][0] = cos(theta2);
-	double theta2 = 0.5 * PI * wvalue/180.0;
+      quat[i][1] = xvalue * sintheta2;
-	double sintheta2 = sin(theta2);
+      quat[i][2] = yvalue * sintheta2;
-	double **quat = atom->quat;
+      quat[i][3] = zvalue * sintheta2;
-	quat[i][0] = cos(theta2);
+      MathExtra::normalize4(quat[i]);
-	quat[i][1] = xvalue * sintheta2;
-	quat[i][2] = yvalue * sintheta2;
-	quat[i][3] = zvalue * sintheta2;
-	MathExtra::normalize4(quat[i]);
-      }
-      count++;
     }
+    count++;
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -397,8 +445,6 @@ void Set::setrandom(int keyword)
 {
   int i;
 
-  if (keyword == DIPOLE_RANDOM) atom->check_dipole();
-
   selection(atom->nlocal);
   RanPark *random = new RanPark(lmp,1);
   double **x = atom->x;
@@ -422,7 +468,6 @@ void Set::setrandom(int keyword)
 
   } else if (keyword == DIPOLE_RANDOM) {
     int *type = atom->type;
-    double *dipole = atom->dipole;
     double **mu = atom->mu;
     int nlocal = atom->nlocal;
 
@@ -431,34 +476,32 @@ void Set::setrandom(int keyword)
     if (domain->dimension == 3) {
       for (i = 0; i < nlocal; i++)
 	if (select[i]) {
-	  if (dipole[type[i]] > 0.0) {
+	  random->reset(seed,x[i]);
-	    random->reset(seed,x[i]);
+	  mu[i][0] = random->uniform() - 0.5;
-	    mu[i][0] = random->uniform() - 0.5;
+	  mu[i][1] = random->uniform() - 0.5;
-	    mu[i][1] = random->uniform() - 0.5;
+	  mu[i][2] = random->uniform() - 0.5;
-	    mu[i][2] = random->uniform() - 0.5;
+	  msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
-	    msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
+	  scale = dvalue/sqrt(msq);
-	    scale = dipole[type[i]]/sqrt(msq);
+	  mu[i][0] *= scale;
-	    mu[i][0] *= scale;
+	  mu[i][1] *= scale;
-	    mu[i][1] *= scale;
+	  mu[i][2] *= scale;
-	    mu[i][2] *= scale;
+	  mu[i][3] = dvalue;
-	    count++;
+	  count++;
-	  }
 	}
 
     } else {
       for (i = 0; i < nlocal; i++)
 	if (select[i]) {
-	  if (dipole[type[i]] > 0.0) {
+	  random->reset(seed,x[i]);
-	    random->reset(seed,x[i]);
+	  mu[i][0] = random->uniform() - 0.5;
-	    mu[i][0] = random->uniform() - 0.5;
+	  mu[i][1] = random->uniform() - 0.5;
-	    mu[i][1] = random->uniform() - 0.5;
+	  mu[i][2] = 0.0;
-	    mu[i][2] = 0.0;
+	  msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1];
-	    msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1];
+	  scale = dvalue/sqrt(msq);
-	    scale = dipole[type[i]]/sqrt(msq);
+	  mu[i][0] *= scale;
-	    mu[i][0] *= scale;
+	  mu[i][1] *= scale;
-	    mu[i][1] *= scale;
+	  mu[i][3] = dvalue;
-	    count++;
+	  count++;
-	  }
 	}
     }
 

src/write_restart.cpp

@@ -56,7 +56,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT,
        SPECIAL_LJ_1,SPECIAL_LJ_2,SPECIAL_LJ_3,
        SPECIAL_COUL_1,SPECIAL_COUL_2,SPECIAL_COUL_3,
        XY,XZ,YZ};
-enum{MASS,SHAPE,DIPOLE};
+enum{MASS};
 enum{PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
 
 enum{IGNORE,WARN,ERROR};                    // same as thermo.cpp
@@ -394,16 +394,6 @@ void WriteRestart::type_arrays()
     fwrite(&flag,sizeof(int),1,fp);
     fwrite(&atom->mass[1],sizeof(double),atom->ntypes,fp);
   }
-  if (atom->shape) {
-    int flag = SHAPE;
-    fwrite(&flag,sizeof(int),1,fp);
-    fwrite(&atom->shape[1][0],sizeof(double),atom->ntypes*3,fp);
-  }
-  if (atom->dipole) {
-    int flag = DIPOLE;
-    fwrite(&flag,sizeof(int),1,fp);
-    fwrite(&atom->dipole[1],sizeof(double),atom->ntypes,fp);
-  }
 
   // -1 flag signals end of type arrays