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

src/ASPHERE/compute_erotate_asphere.cpp

@@ -15,6 +15,7 @@
 #include "compute_erotate_asphere.h"
 #include "math_extra.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "update.h"
 #include "force.h"
 #include "memory.h"
@@ -24,19 +25,27 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-ComputeERotateAsphere::ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) :
+ComputeERotateAsphere::
+ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
   if (narg != 3) error->all("Illegal compute erotate/asphere command");
 
-  if (!atom->angmom_flag || !atom->quat_flag)
-    error->all("Compute erotate/asphere requires atom attributes angmom, quat");
-
   scalar_flag = 1;
   extscalar = 1;
 
   inertia = 
     memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia");
+
+  // error checks
+
+  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");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -50,11 +59,21 @@ ComputeERotateAsphere::~ComputeERotateAsphere()
 
 void ComputeERotateAsphere::init()
 {
+  // check that all particles are finite-size
+  // no point particles allowed, spherical is OK
+
+  double **shape = atom->shape;
+  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)
+      if (shape[type[i]][0] == 0.0)
+	error->one("Compue erotate/asphere requires extended particles");
+
   pfactor = 0.5 * force->mvv2e;
-
-  if (!atom->shape)
-    error->all("Compute erotate/asphere requires atom attribute shape");
-
   calculate_inertia();
 }
 
@@ -62,6 +81,8 @@ void ComputeERotateAsphere::init()
 
 double ComputeERotateAsphere::compute_scalar()
 {
+  int i,itype;
+
   invoked_scalar = update->ntimestep;
 
   double **quat = atom->quat;
@@ -70,12 +91,14 @@ double ComputeERotateAsphere::compute_scalar()
   int *type = atom->type;
   int nlocal = atom->nlocal;
 
-  int itype;
+  // sum rotational energy for each particle
+  // no point particles since divide by inertia
+
   double wbody[3];
   double rot[3][3];
   double erotate = 0.0;
 
-  for (int i = 0; i < nlocal; i++)
+  for (i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       itype = type[i];
 

src/ASPHERE/compute_temp_asphere.cpp

@@ -19,6 +19,7 @@
 #include "compute_temp_asphere.h"
 #include "math_extra.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "update.h"
 #include "force.h"
 #include "domain.h"
@@ -38,9 +39,6 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
   if (narg != 3 && narg != 4)
     error->all("Illegal compute temp/asphere command");
 
-  if (!atom->quat_flag || !atom->angmom_flag)
-    error->all("Compute temp/asphere requires atom attributes quat, angmom");
-
   scalar_flag = vector_flag = 1;
   size_vector = 6;
   extscalar = 0;
@@ -59,6 +57,16 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
   vector = new double[6];
   inertia = 
     memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia");
+
+  // error checks
+
+  if (!atom->angmom_flag || !atom->quat_flag ||
+      !atom->avec->shape_type)
+    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,6 +82,20 @@ ComputeTempAsphere::~ComputeTempAsphere()
 
 void ComputeTempAsphere::init()
 {
+  // check that all particles are finite-size
+  // no point particles allowed, spherical is OK
+
+  double **shape = atom->shape;
+  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)
+      if (shape[type[i]][0] == 0.0)
+	error->one("Compue temp/asphere requires extended particles");
+
   if (tempbias) {
     int i = modify->find_compute(id_bias);
     if (i < 0) error->all("Could not find compute ID for temperature bias");
@@ -101,44 +123,30 @@ void ComputeTempAsphere::init()
 
 void ComputeTempAsphere::dof_compute()
 {
+  // 6 dof for 3d, 3 dof for 2d
+  // assume full rotation of extended particles
+  // user can correct this via compute_modify if needed
+
   double natoms = group->count(igroup);
-  int dimension = domain->dimension;
+  if (domain->dimension == 3) dof = 6*natoms;
-  dof = dimension * natoms;
+  else dof = 3*natoms;
+
+  // additional adjustments to dof
 
   if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms;
+  else if (tempbias == 2) {
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+    int count = 0;
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	if (tbias->dof_remove(i)) count++;
+    int count_all;
+    MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world);
+    if (domain->dimension == 3) dof -= 6*count_all;
+    else dof -= 3*count_all;
+  }
 
-  // rotational degrees of freedom
-  // 0 for sphere, 2 for uniaxial, 3 for biaxial
-
-  double **shape = atom->shape;
-  int *type = atom->type;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  int itype;
-  int count = 0;
-
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) {
-      if (tempbias == 2 && tbias->dof_remove(i)) continue;
-      itype = type[i];
-      if (dimension == 2) {
-	if (shape[itype][0] == shape[itype][1]) continue;
-	else count++;
-      } else {
-	if (shape[itype][0] == shape[itype][1] && 
-	    shape[itype][1] == shape[itype][2]) continue;
-	else if (shape[itype][0] == shape[itype][1] || 
-		 shape[itype][1] == shape[itype][2] ||
-		 shape[itype][0] == shape[itype][2]) count += 2;
-	else count += 3;
-      }
-    }
-
-  int count_all;
-  MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world);
-  dof += count_all;
-    
   dof -= extra_dof + fix_dof;
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
@@ -169,31 +177,26 @@ double ComputeTempAsphere::compute_scalar()
   double rot[3][3];
   double t = 0.0;
   
+  // sum translationals and rotational energy for each particle
+  // no point particles since divide by inertia
+
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
 
-      // translational kinetic energy
-      
       itype = type[i];
       t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[itype];
 
       // 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[1] /= inertia[itype][1];
+      wbody[2] /= inertia[itype][2];
 
-      if (!(shape[itype][0] == shape[itype][1] && 
+      t += inertia[itype][0]*wbody[0]*wbody[0] +
-            shape[itype][1] == shape[itype][2])) {
+	inertia[itype][1]*wbody[1]*wbody[1] +
-
+	inertia[itype][2]*wbody[2]*wbody[2];
-        MathExtra::quat_to_mat(quat[i],rot);
-        MathExtra::transpose_times_column3(rot,angmom[i],wbody);
-        wbody[0] /= inertia[itype][0];
-        wbody[1] /= inertia[itype][1];
-        wbody[2] /= inertia[itype][2];
-
-        // rotational kinetic energy
-
-        t += inertia[itype][0]*wbody[0]*wbody[0]+
-             inertia[itype][1]*wbody[1]*wbody[1]+
-             inertia[itype][2]*wbody[2]*wbody[2];
-      }
     }
 
   if (tempbias) tbias->restore_bias_all();

src/ASPHERE/fix_npt_asphere.cpp

@@ -27,6 +27,7 @@
 #include "kspace.h"
 #include "update.h"
 #include "domain.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -38,23 +39,56 @@ enum{NOBIAS,BIAS};
 FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) :
   FixNPT(lmp, narg, arg)
 {
+  inertia = 
+    memory->create_2d_double_array(atom->ntypes+1,3,"fix_npt_asphere:inertia");
+
+  // error checks
+
   if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag ||
-      !atom->shape)
+      !atom->avec->shape_type)
     error->all("Fix npt/asphere requires atom attributes "
 	       "quat, angmom, torque, shape");
+  if (atom->radius_flag || atom->rmass_flag)
+    error->all("Fix npt/asphere cannot be used with atom attributes "
+	       "diameter or rmass");
 }
 
-/* ----------------------------------------------------------------------
+/* ---------------------------------------------------------------------- */
-   1st half of Verlet update 
+
-------------------------------------------------------------------------- */
+FixNPTAsphere::~FixNPTAsphere()
+{
+  memory->destroy_2d_double_array(inertia);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNPTAsphere::init()
+{
+  // check that all particles are finite-size
+  // no point particles allowed, spherical is OK
+
+  double **shape = atom->shape;
+  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)
+      if (shape[type[i]][0] == 0.0)
+	error->one("Fix nvt/asphere requires extended particles");
+
+  FixNPT::init();
+  calculate_inertia();
+}
+
+/* ---------------------------------------------------------------------- */
 
 void FixNPTAsphere::initial_integrate(int vflag)
 {
   int i;
   double dtfm;
 
-  dtq = 0.5 * dtv;
-
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
 
@@ -85,8 +119,6 @@ void FixNPTAsphere::initial_integrate(int vflag)
   }
   factor_rotate = exp(-dthalf*eta_dot);
 
-  // update v of only atoms in group
-
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
@@ -135,7 +167,11 @@ void FixNPTAsphere::initial_integrate(int vflag)
     }
   }
 
-  // update angular momentum by 1/2 step
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dtq = 0.5 * dtv;
+
+  // update angular momentum by 1/2 step for all particles
   // update quaternion a full step via Richardson iteration
   // returns new normalized quaternion
 
@@ -145,9 +181,7 @@ void FixNPTAsphere::initial_integrate(int vflag)
       angmom[i][1] = angmom[i][1]*factor_rotate + dtf*torque[i][1];
       angmom[i][2] = angmom[i][2]*factor_rotate + dtf*torque[i][2];
 		
-      double inertia[3];
+      richardson(quat[i],angmom[i],inertia[type[i]]);
-      calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
-      richardson(quat[i],angmom[i],inertia);
     }
   }
 
@@ -158,9 +192,7 @@ void FixNPTAsphere::initial_integrate(int vflag)
   if (kspace_flag) force->kspace->setup();
 }
 
-/* ----------------------------------------------------------------------
+/* ---------------------------------------------------------------------- */
-   2nd half of Verlet update 
-------------------------------------------------------------------------- */
 
 void FixNPTAsphere::final_integrate()
 {
@@ -320,15 +352,22 @@ void FixNPTAsphere::omega_from_mq(double *q, double *m, double *inertia,
   MathExtra::times_column3(rot,wbody,w);
 }
 
+
 /* ----------------------------------------------------------------------
-   calculate the moment of inertia for an ellipsoid, from mass and radii
+   principal moments of inertia for ellipsoids
-   shape = x,y,z radii in body frame
 ------------------------------------------------------------------------- */
 
-void FixNPTAsphere::calculate_inertia(double mass, double *shape,
+void FixNPTAsphere::calculate_inertia()
-				      double *inertia)
 {
-  inertia[0] = 0.2*mass*(shape[1]*shape[1]+shape[2]*shape[2]);
+  double *mass = atom->mass;
-  inertia[1] = 0.2*mass*(shape[0]*shape[0]+shape[2]*shape[2]);
+  double **shape = atom->shape;
-  inertia[2] = 0.2*mass*(shape[0]*shape[0]+shape[1]*shape[1]);
+  
+  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_npt_asphere.h

@@ -21,16 +21,19 @@ namespace LAMMPS_NS {
 class FixNPTAsphere : public FixNPT {
  public:
   FixNPTAsphere(class LAMMPS *, int, char **);
+  ~FixNPTAsphere();
+  void init();
   void initial_integrate(int);
   void final_integrate();
 
  private:
   double dtq;
   double factor_rotate;
+  double **inertia;
 
   void richardson(double *, double *, double *);
   void omega_from_mq(double *, double *, double *, double *);
-  void calculate_inertia(double mass, double *shape, double *inertia);
+  void calculate_inertia();
 };
 
 }

src/ASPHERE/fix_nve_asphere.cpp

@@ -37,12 +37,18 @@ using namespace LAMMPS_NS;
 FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) : 
   FixNVE(lmp, narg, arg)
 {
-  if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag ||
-      !atom->shape)
-    error->all("Fix nve/asphere requires atom attributes "
-	       "quat, angmom, torque, shape");
   inertia = 
-    memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia");
+    memory->create_2d_double_array(atom->ntypes+1,3,"fix_nve_asphere:inertia");
+
+  // error checks
+
+  if (!atom->angmom_flag || !atom->quat_flag || !atom->torque_flag ||
+      !atom->avec->shape_type)
+    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");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -56,6 +62,20 @@ FixNVEAsphere::~FixNVEAsphere()
 
 void FixNVEAsphere::init()
 {
+  // check that all particles are finite-size
+  // no point particles allowed, spherical is OK
+
+  double **shape = atom->shape;
+  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)
+      if (shape[type[i]][0] == 0.0)
+	error->one("Fix nve/asphere requires extended particles");
+
   FixNVE::init();
   calculate_inertia();
 }
@@ -66,8 +86,6 @@ void FixNVEAsphere::initial_integrate(int vflag)
 {
   double dtfm;
 
-  dtq = 0.5 * dtv;
-
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
@@ -80,6 +98,10 @@ void FixNVEAsphere::initial_integrate(int vflag)
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dtq = 0.5 * dtv;
+
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       dtfm = dtf / mass[type[i]];
@@ -216,7 +238,7 @@ 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]);

src/ASPHERE/fix_nvt_asphere.cpp

@@ -28,6 +28,7 @@
 #include "update.h"
 #include "modify.h"
 #include "compute.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -39,10 +40,47 @@ enum{NOBIAS,BIAS};
 FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) :
   FixNVT(lmp, narg, arg)
 {
+  inertia = 
+    memory->create_2d_double_array(atom->ntypes+1,3,"fix_nvt_asphere:inertia");
+
+  // error checks
+
   if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag ||
-      !atom->shape)
+      !atom->avec->shape_type)
     error->all("Fix nvt/asphere requires atom attributes "
 	       "quat, angmom, torque, shape");
+  if (atom->radius_flag || atom->rmass_flag)
+    error->all("Fix nvt/asphere cannot be used with atom attributes "
+	       "diameter or rmass");
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixNVTAsphere::~FixNVTAsphere()
+{
+  memory->destroy_2d_double_array(inertia);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTAsphere::init()
+{
+  // check that all particles are finite-size
+  // no point particles allowed, spherical is OK
+
+  double **shape = atom->shape;
+  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)
+      if (shape[type[i]][0] == 0.0)
+	error->one("Fix nvt/asphere requires extended particles");
+
+  FixNVT::init();
+  calculate_inertia();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -51,8 +89,6 @@ void FixNVTAsphere::initial_integrate(int vflag)
 {
   double dtfm;
 
-  dtq = 0.5 * dtv;
-
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
   t_target = t_start + delta * (t_stop-t_start);
@@ -65,8 +101,6 @@ void FixNVTAsphere::initial_integrate(int vflag)
   eta += dtv*eta_dot;
   factor = exp(-dthalf*eta_dot);
 
-  // update v and x of only atoms in group
-
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
@@ -79,6 +113,10 @@ void FixNVTAsphere::initial_integrate(int vflag)
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dtq = 0.5 * dtv;
+
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -98,9 +136,7 @@ void FixNVTAsphere::initial_integrate(int vflag)
 	angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1];
 	angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2];
 	
-	double inertia[3];
+	richardson(quat[i],angmom[i],inertia[type[i]]);
-	calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
-	richardson(quat[i],angmom[i],inertia);
       }
     }
 
@@ -125,9 +161,7 @@ void FixNVTAsphere::initial_integrate(int vflag)
 	angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1];
 	angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2];
 	
-	double inertia[3];
+	richardson(quat[i],angmom[i],inertia[type[i]]);
-	calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
-	richardson(quat[i],angmom[i],inertia);
       }
     }
   }
@@ -139,8 +173,6 @@ void FixNVTAsphere::final_integrate()
 {
   double dtfm;
 
-  // update v of only atoms in group
-
   double **v = atom->v;
   double **f = atom->f;
   double **angmom = atom->angmom;
@@ -269,14 +301,20 @@ void FixNVTAsphere::omega_from_mq(double *q, double *m, double *inertia,
 }
 
 /* ----------------------------------------------------------------------
-   calculate the moment of inertia for an ELLIPSOID, from mass and radii
+   principal moments of inertia for ellipsoids
-   shape = x,y,z radii in body frame
 ------------------------------------------------------------------------- */
 
-void FixNVTAsphere::calculate_inertia(double mass, double *shape,
+void FixNVTAsphere::calculate_inertia()
-				      double *inertia)
 {
-  inertia[0] = 0.2*mass*(shape[1]*shape[1]+shape[2]*shape[2]);
+  double *mass = atom->mass;
-  inertia[1] = 0.2*mass*(shape[0]*shape[0]+shape[2]*shape[2]);
+  double **shape = atom->shape;
-  inertia[2] = 0.2*mass*(shape[0]*shape[0]+shape[1]*shape[1]);
+  
+  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_nvt_asphere.h

@@ -21,15 +21,18 @@ namespace LAMMPS_NS {
 class FixNVTAsphere : public FixNVT {
  public:
   FixNVTAsphere(class LAMMPS *, int, char **);
+  ~FixNVTAsphere();
+  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(double mass, double *shape, double *inertia);
+  void calculate_inertia();
 };
 
 }

src/ASPHERE/pair_gayberne.cpp

@@ -22,6 +22,7 @@
 #include "pair_gayberne.h"
 #include "math_extra.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
@@ -334,8 +335,10 @@ void PairGayBerne::coeff(int narg, char **arg)
 
 void PairGayBerne::init_style()
 {
-  if (!atom->quat_flag || !atom->torque_flag)
+  if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type)
-    error->all("Pair gayberne requires atom attributes quat, torque");
+    error->all("Pair gayberne requires atom attributes quat, torque, shape");
+  if (atom->radius_flag)
+    error->all("Pair gayberne cannot be used with atom attribute diameter");
 
   int irequest = neighbor->request(this);
 

src/ASPHERE/pair_resquared.cpp

@@ -22,6 +22,7 @@
 #include "pair_resquared.h"
 #include "math_extra.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
@@ -323,8 +324,10 @@ void PairRESquared::coeff(int narg, char **arg)
 
 void PairRESquared::init_style()
 {
-  if (!atom->quat_flag || !atom->torque_flag)
+  if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type)
-    error->all("Pair resquared requires atom attributes quat, torque");
+    error->all("Pair resquared requires atom attributes quat, torque, shape");
+  if (atom->radius_flag)
+    error->all("Pair resquared cannot be used with atom attribute diameter");
 
   int irequest = neighbor->request(this);
 
@@ -353,9 +356,11 @@ double PairRESquared::init_one(int i, int j)
     error->all("Pair resquared epsilon a,b,c coeffs are not all set");
 
   int ishape = 0;
-  if (shape[i][0] != 0 && shape[i][1] != 0 && shape[i][2] != 0) ishape = 1;
+  if (shape[i][0] != 0.0 && shape[i][1] != 0.0 && shape[i][2] != 0.0) 
+    ishape = 1;
   int jshape = 0;
-  if (shape[j][0] != 0 && shape[j][1] != 0 && shape[j][2] != 0) jshape = 1;
+  if (shape[j][0] != 0.0 && shape[j][1] != 0.0 && shape[j][2] != 0.0) 
+    jshape = 1;
   
   if (ishape == 0 && jshape == 0) {
     form[i][j] = SPHERE_SPHERE;

src/COLLOID/pair_lubricate.cpp

@@ -21,6 +21,7 @@
 #include "string.h"
 #include "pair_lubricate.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
@@ -56,6 +57,7 @@ PairLubricate::~PairLubricate()
     memory->destroy_2d_double_array(cut);
     memory->destroy_2d_double_array(cut_inner);
   }
+
   delete random;
 }
 
@@ -99,10 +101,10 @@ void PairLubricate::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  a_squeeze = a_shear = a_pump = a_twist = 0.0;
-
   // loop over neighbors of my atoms
 
+  a_squeeze = a_shear = a_pump = a_twist = 0.0;
+
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     xtmp = x[i][0];
@@ -388,17 +390,21 @@ void PairLubricate::coeff(int narg, char **arg)
 
 void PairLubricate::init_style()
 {
-  if (!atom->torque_flag || atom->shape == NULL)
+  if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type)
-    error->all("Pair lubricate requires atom attributes torque, shape");
+    error->all("Pair lubricate requires atom attributes quat, torque, shape");
+  if (atom->radius_flag || atom->rmass_flag)
+    error->all("Pair lubricate cannot be used with atom attributes"
+	       "diameter or rmass");
 
-  // insure all particle shapes are spherical
+  // insure all particle shapes are finite-size, spherical, and monodisperse
 
+  double value = atom->shape[1][0];
+  if (value == 0.0) error->all("Pair lubricate requires extended particles");
   for (int i = 1; i <= atom->ntypes; i++)
-    if ((atom->shape[i][0] != atom->shape[i][1]) ||
+    if (atom->shape[i][0] != value || atom->shape[i][0] != value || 
-	(atom->shape[i][0] != atom->shape[i][2]) ||
+	atom->shape[i][0] != value)
-	(atom->shape[i][1] != atom->shape[i][2]) )
+      error->all("Pair lubricate requires spherical, mono-disperse particles");
-      error->all("Pair lubricate requires spherical particles");
+  
-
   int irequest = neighbor->request(this);
 }
 

src/GRANULAR/atom_vec_granular.cpp

@@ -533,8 +533,11 @@ int AtomVecGranular::unpack_restart(double *buf)
 
   radius[nlocal] = buf[m++];
   density[nlocal] = buf[m++];
-  rmass[nlocal] = 4.0*PI/3.0 * 
+
-    radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
+  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++];
@@ -601,10 +604,17 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values)
     error->one("Invalid atom type in Atoms section of data file");
 
   radius[nlocal] = 0.5 * atof(values[2]);
+  if (radius[nlocal] < 0.0)
+    error->one("Invalid radius in Atoms section of data file");
+
   density[nlocal] = atof(values[3]);
-  rmass[nlocal] = 4.0*PI/3.0 *
+  if (density[nlocal] <= 0.0)
-    radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
+    error->one("Invalid density in Atoms section of data file");
-  if (rmass[nlocal] <= 0.0) error->one("Invalid mass value");
+
+  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];
 
   x[nlocal][0] = coord[0];
   x[nlocal][1] = coord[1];
@@ -631,10 +641,17 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values)
 int AtomVecGranular::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]);
-  rmass[nlocal] = 4.0*PI/3.0 *
+  if (density[nlocal] <= 0.0)
-    radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
+    error->one("Invalid density in Atoms section of data file");
-  if (rmass[nlocal] <= 0.0) error->one("Invalid mass value");
+
+  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];
 
   v[nlocal][0] = 0.0;
   v[nlocal][1] = 0.0;

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -43,7 +43,7 @@ PairGranHertzHistory::PairGranHertzHistory(LAMMPS *lmp) :
 
 void PairGranHertzHistory::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum;
+  int i,j,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
   double radi,radj,radsum,rsq,r,rinv,rsqinv;
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -66,6 +66,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
   double **torque = atom->torque;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
@@ -141,9 +143,18 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
 
 	// normal force = Hertzian contact + normal velocity damping
 
-	meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+	if (rmass) {
-	if (mask[i] & freeze_group_bit) meff = rmass[j];
+	  meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
-	if (mask[j] & freeze_group_bit) meff = rmass[i];
+	  if (mask[i] & freeze_group_bit) meff = rmass[j];
+	  if (mask[j] & freeze_group_bit) meff = rmass[i];
+	} else {
+	  itype = type[i];
+	  jtype = type[j];
+	  meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+	  if (mask[i] & freeze_group_bit) meff = mass[jtype];
+	  if (mask[j] & freeze_group_bit) meff = mass[itype];
+	}
+
 	damp = meff*gamman*vnnr*rsqinv;
 	ccel = kn*(radsum-r)*rinv - damp;
 	polyhertz = sqrt((radsum-r)*radi*radj / radsum);

src/GRANULAR/pair_gran_hooke.cpp

@@ -39,7 +39,7 @@ PairGranHooke::PairGranHooke(LAMMPS *lmp) : PairGranHookeHistory(lmp)
 
 void PairGranHooke::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum;
+  int i,j,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
   double radi,radj,radsum,rsq,r,rinv,rsqinv;
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -59,6 +59,8 @@ void PairGranHooke::compute(int eflag, int vflag)
   double **torque = atom->torque;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
@@ -121,9 +123,18 @@ void PairGranHooke::compute(int eflag, int vflag)
 
 	// normal forces = Hookian contact + normal velocity damping
 
-	meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+	if (rmass) {
-	if (mask[i] & freeze_group_bit) meff = rmass[j];
+	  meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
-	if (mask[j] & freeze_group_bit) meff = rmass[i];
+	  if (mask[i] & freeze_group_bit) meff = rmass[j];
+	  if (mask[j] & freeze_group_bit) meff = rmass[i];
+	} else {
+	  itype = type[i];
+	  jtype = type[j];
+	  meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+	  if (mask[i] & freeze_group_bit) meff = mass[jtype];
+	  if (mask[j] & freeze_group_bit) meff = mass[itype];
+	}
+
 	damp = meff*gamman*vnnr*rsqinv;
 	ccel = kn*(radsum-r)*rinv - damp;
 

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -67,7 +67,7 @@ PairGranHookeHistory::~PairGranHookeHistory()
 
 void PairGranHookeHistory::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum;
+  int i,j,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
   double radi,radj,radsum,rsq,r,rinv,rsqinv;
   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -90,6 +90,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
   double **torque = atom->torque;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
@@ -165,9 +167,18 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 
 	// normal forces = Hookian contact + normal velocity damping
 
-	meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+	if (rmass) {
-	if (mask[i] & freeze_group_bit) meff = rmass[j];
+	  meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
-	if (mask[j] & freeze_group_bit) meff = rmass[i];
+	  if (mask[i] & freeze_group_bit) meff = rmass[j];
+	  if (mask[j] & freeze_group_bit) meff = rmass[i];
+	} else {
+	  itype = type[i];
+	  jtype = type[j];
+	  meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+	  if (mask[i] & freeze_group_bit) meff = mass[jtype];
+	  if (mask[j] & freeze_group_bit) meff = mass[itype];
+	}
+
 	damp = meff*gamman*vnnr*rsqinv;
 	ccel = kn*(radsum-r)*rinv - damp;
 
@@ -328,6 +339,8 @@ void PairGranHookeHistory::init_style()
 {
   int i;
 
+  // error and warning checks
+
   if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)
     error->all("Pair granular requires atom attributes radius, omega, torque");
   if (atom->avec->ghost_velocity == 0)

src/atom.cpp

@@ -1124,17 +1124,15 @@ void Atom::set_shape(const char *str)
   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)
+  if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || 
-    error->all("Invalid shape value");
+      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");
-  if (shape[itype][1] == 0.0 &&
-      (shape[itype][0] != 0.0 || shape[itype][2] != 0.0))
-    error->all("Invalid shape value");
-  if (shape[itype][2] == 0.0 &&
-      (shape[itype][0] != 0.0 || shape[itype][1] != 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");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1159,18 +1157,15 @@ void Atom::set_shape(int narg, char **arg)
     shape[itype][2] = 0.5*atof(arg[3]);
     shape_setflag[itype] = 1;
 
-    if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 ||
+    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 &&
+    if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || 
-	(shape[itype][1] != 0.0 || shape[itype][2] != 0.0))
+	shape[itype][2] > 0.0) {
-      error->all("Invalid shape value");
+      if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || 
-    if (shape[itype][1] == 0.0 &&
+	  shape[itype][2] == 0.0)
-	(shape[itype][0] != 0.0 || shape[itype][2] != 0.0))
+	error->all("Invalid shape value");
-      error->all("Invalid shape value");
+    }
-    if (shape[itype][2] == 0.0 &&
-	(shape[itype][0] != 0.0 || shape[itype][1] != 0.0))
-      error->all("Invalid shape value");
   }
 }
 

src/compute_erotate_sphere.cpp

@@ -14,6 +14,7 @@
 #include "mpi.h"
 #include "compute_erotate_sphere.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "update.h"
 #include "force.h"
 #include "domain.h"
@@ -31,71 +32,100 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 3) error->all("Illegal compute erotate/sphere command");
 
-  if (!atom->omega_flag) 
-    error->all("Compute erotate/sphere requires atom attribute omega");
-
   scalar_flag = 1;
   extscalar = 1;
 
-  inertia = new double[atom->ntypes+1];
+  // error checks
-}
 
-/* ---------------------------------------------------------------------- */
+  if (!atom->omega_flag) 
-
+    error->all("Compute erotate/sphere requires atom attribute omega");
-ComputeERotateSphere::~ComputeERotateSphere()
+  if (!atom->radius_flag && !atom->avec->shape_type)
-{
+    error->all("Compute erotate/sphere requires atom attribute "
-  delete [] inertia;
+	       "radius or shape");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeERotateSphere::init()
 {
-  pfactor = 0.5 * force->mvv2e * INERTIA;
+  int i,itype;
 
-  if (atom->mass && !atom->shape)
+  // if shape used, check that all particles are spherical
-    error->all("Compute erotate/sphere requires atom attribute shape");
+  // point particles are allowed
-  if (atom->rmass && !atom->radius_flag)
-    error->all("Compute erotate/sphere requires atom attribute radius");
 
-  if (atom->mass) {
+  if (atom->radius == NULL) {
-    double *mass = atom->mass;
     double **shape = atom->shape;
-    
+    int *type = atom->type;
-    for (int i = 1; i <= atom->ntypes; i++) {
+    int *mask = atom->mask;
-      if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
+    int nlocal = atom->nlocal;
-	error->all("Compute erotate/sphere requires spherical particle shapes");
+
-      inertia[i] = shape[i][0]*shape[i][0] * mass[i];
+    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;
 }
 
 /* ---------------------------------------------------------------------- */
 
 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
+
   double erotate = 0.0;
 
-  if (rmass) {
+  if (radius) {
-    for (int i = 0; i < nlocal; i++) 
+    if (rmass) {
-      if (mask[i] & groupbit)
+      for (i = 0; i < nlocal; i++) 
-	erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
+	if (mask[i] & groupbit)
-		    omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];
+	  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)
+	  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 {
-    for (int i = 0; i < nlocal; i++)
+    if (rmass) {
-      if (mask[i] & groupbit)
+      for (i = 0; i < nlocal; i++) 
-	erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
+	if (mask[i] & groupbit) {
-		    omega[i][2]*omega[i][2]) * inertia[type[i]];
+	  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);

src/compute_erotate_sphere.h

@@ -21,13 +21,12 @@ namespace LAMMPS_NS {
 class ComputeERotateSphere : public Compute {
  public:
   ComputeERotateSphere(class LAMMPS *, int, char **);
-  ~ComputeERotateSphere();
+  ~ComputeERotateSphere() {}
   void init();
   double compute_scalar();
 
  private:
   double pfactor;
-  double *inertia;
 };
 
 }

src/compute_group_group.cpp

@@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing authors: Naveen Michaud-Agrawal (Johns Hopkins U)
+   Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
 ------------------------------------------------------------------------- */
 
 #include "mpi.h"

src/compute_temp_sphere.cpp

@@ -15,6 +15,7 @@
 #include "string.h"
 #include "compute_temp_sphere.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "update.h"
 #include "force.h"
 #include "domain.h"
@@ -35,9 +36,6 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
   if (narg != 3 && narg != 4)
     error->all("Illegal compute temp/sphere command");
 
-  if (!atom->omega_flag)
-    error->all("Compute temp/sphere requires atom attribute omega");
-
   scalar_flag = vector_flag = 1;
   size_vector = 6;
   extscalar = 0;
@@ -54,23 +52,51 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
   }
 
   vector = new double[6];
-  inertia = new double[atom->ntypes+1];
+
+  // error checks
+
+  if (!atom->omega_flag)
+    error->all("Compute temp/sphere requires atom attribute omega");
+  if (!atom->radius_flag && !atom->avec->shape_type)
+    error->all("Compute temp/sphere requires atom attribute "
+	       "radius or shape");
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputeTempSphere::~ComputeTempSphere()
 {
+  delete [] id_bias;
   delete [] vector;
-  delete [] inertia;
 }
 
 /* ---------------------------------------------------------------------- */
 
 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) {
-    int i = modify->find_compute(id_bias);
+    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)
@@ -85,44 +111,119 @@ void ComputeTempSphere::init()
   }
 
   fix_dof = 0;
-  for (int i = 0; i < modify->nfix; i++)
+  for (i = 0; i < modify->nfix; i++)
     fix_dof += modify->fix[i]->dof(igroup);
   dof_compute();
-
-  if (atom->mass) {
-    double *mass = atom->mass;
-    double **shape = atom->shape;
-    
-    for (int i = 1; i <= atom->ntypes; i++) {
-      if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
-	error->all("Compute temp/sphere requires spherical particle shapes");
-      inertia[i] = INERTIA * shape[i][0]*shape[i][0] * mass[i];
-    }
-  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeTempSphere::dof_compute()
 {
-  double natoms = group->count(igroup);
+  int count,count_all;
-  int nper = 6;
-  if (domain->dimension == 2) nper = 3;
-  dof = nper * natoms;
 
-  if (tempbias) {
+  // 6 or 3 dof for extended/point particles for 3d
-    if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms;
+  // 3 or 2 dof for extended/point particles for 2d
-    else {
+  // assume full rotation of extended particles
-      int *mask = atom->mask;
+  // user can correct this via compute_modify if needed
-      int nlocal = atom->nlocal;
+
-      int count = 0;
+  int dimension = domain->dimension;
-      for (int i = 0; i < nlocal; i++)
+
-	if (mask[i] & groupbit)
+  double *radius = atom->radius;
-	  if (tbias->dof_remove(i)) count++;
+  double **shape = atom->shape;
-      int count_all;
+  int *type = atom->type;
-      MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world);
+  int *mask = atom->mask;
-      dof -= nper * count_all;
+  int nlocal = atom->nlocal;
+
+  count = 0;
+  if (dimension == 3) {
+    if (radius) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  if (radius[i] == 0.0) count += 3;
+	  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++) {
+	if (mask[i] & groupbit) {
+	  if (radius[i] == 0.0) count += 2;
+	  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);
+  dof = count_all;
+
+  // additional adjustments to dof
+
+  if (tempbias == 1) {
+    double natoms = group->count(igroup);
+    dof -= tbias->dof_remove(-1) * natoms;
+
+  } else if (tempbias == 2) {
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    count = 0;
+    if (dimension == 3) {
+      if (radius) {
+	for (int i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit) {
+	    if (tbias->dof_remove(i)) {
+	      if (radius[i] == 0.0) count += 3;
+	      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++)
+	  if (mask[i] & groupbit) {
+	    if (tbias->dof_remove(i)) {
+	      if (radius[i] == 0.0) count += 2;
+	      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);
+    dof -= count_all;
   }
 
   dof -= extra_dof + fix_dof;
@@ -134,6 +235,8 @@ void ComputeTempSphere::dof_compute()
 
 double ComputeTempSphere::compute_scalar()
 {
+  int i,itype;
+
   invoked_scalar = update->ntimestep;
 
   if (tempbias) {
@@ -143,30 +246,65 @@ double ComputeTempSphere::compute_scalar()
 
   double **v = atom->v;
   double **omega = atom->omega;
-  double *mass = atom->mass;
-  double *rmass = atom->rmass;
   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 or shape = 0
+
   double t = 0.0;
 
-  if (rmass) {
+  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 {
-    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]) * 
+	if (mask[i] & groupbit) {
-	  mass[type[i]];
+	  itype = type[i];
-	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[type[i]];
+	    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();
@@ -181,7 +319,7 @@ double ComputeTempSphere::compute_scalar()
 
 void ComputeTempSphere::compute_vector()
 {
-  int i;
+  int i,itype;
 
   invoked_vector = update->ntimestep;
 
@@ -195,51 +333,103 @@ void ComputeTempSphere::compute_vector()
   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 or shape = 0
+
   double massone,inertiaone,t[6];
   for (i = 0; i < 6; i++) t[i] = 0.0;
 
-  if (rmass) {
+  if (radius) {
-    for (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];
+	  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*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];
+	}
+    }
 
-	inertiaone = INERTIA*radius[i]*radius[i]*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 (rmass) {
-      if (mask[i] & groupbit) {
+      for (i = 0; i < nlocal; i++)
-	massone = mass[type[i]];
+	if (mask[i] & groupbit) {
-	t[0] += massone * v[i][0]*v[i][0];
+	  itype = type[i];
-	t[1] += massone * v[i][1]*v[i][1];
+	  massone = rmass[i];
-	t[2] += massone * v[i][2]*v[i][2];
+	  t[0] += massone * v[i][0]*v[i][0];
-	t[3] += massone * v[i][0]*v[i][1];
+	  t[1] += massone * v[i][1]*v[i][1];
-	t[4] += massone * v[i][0]*v[i][2];
+	  t[2] += massone * v[i][2]*v[i][2];
-	t[5] += massone * v[i][1]*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];
+	}
 
-	inertiaone = inertia[type[i]];
+    } else {
-	t[0] += inertiaone * omega[i][0]*omega[i][0];
+      for (i = 0; i < nlocal; i++)
-	t[1] += inertiaone * omega[i][1]*omega[i][1];
+	if (mask[i] & groupbit) {
-	t[2] += inertiaone * omega[i][2]*omega[i][2];
+	  itype = type[i];
-	t[3] += inertiaone * omega[i][0]*omega[i][1];
+	  massone = mass[itype];
-	t[4] += inertiaone * omega[i][0]*omega[i][2];
+	  t[0] += massone * v[i][0]*v[i][0];
-	t[5] += inertiaone * omega[i][1]*omega[i][2];
+	  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();

src/fix_nph.cpp

@@ -250,8 +250,6 @@ int FixNPH::setmask()
 void FixNPH::init()
 {
   if (domain->triclinic) error->all("Cannot use fix nph with triclinic box");
-  if (atom->mass == NULL)
-    error->all("Cannot use fix nph without per-type mass defined");
 
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"deform") == 0) {
@@ -353,6 +351,7 @@ void FixNPH::setup(int vflag)
 void FixNPH::initial_integrate(int vflag)
 {
   int i;
+  double dtfm;
 
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
@@ -380,18 +379,29 @@ void FixNPH::initial_integrate(int vflag)
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double dtfm;
+  if (rmass) {
-  for (i = 0; i < nlocal; i++) {
+    for (i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
+      if (mask[i] & groupbit) {
-      dtfm = dtf / mass[type[i]];
+	dtfm = dtf / rmass[i];
-      v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-      v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
-      v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+      }
+    }
+  } else {
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	dtfm = dtf / mass[type[i]];
+	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+      }
     }
   }
 
@@ -423,23 +433,35 @@ void FixNPH::initial_integrate(int vflag)
 void FixNPH::final_integrate()
 {
   int i;
+  double dtfm;
 
   // v update only for atoms in NPH group
 
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double dtfm;
+  if (rmass) {
-  for (i = 0; i < nlocal; i++) {
+    for (i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
+      if (mask[i] & groupbit) {
-      dtfm = dtf / mass[type[i]];
+	dtfm = dtf / rmass[i];
-      v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
-      v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
-      v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+      }
+    }
+  } else {
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	dtfm = dtf / mass[type[i]];
+	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+      }
     }
   }
 
@@ -499,6 +521,7 @@ void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag)
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -533,12 +556,23 @@ void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag)
 
     // v update only for atoms in NPH group
 
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  dtfm = dtf / rmass[i];
-	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
+      }
+    } else {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
       }
     }
 
@@ -550,12 +584,23 @@ void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag)
 
     // v update only for atoms in NPH group
 
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[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];
+	}
+      }
+    } else {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] += dtfm*f[i][0];
+	  v[i][1] += dtfm*f[i][1];
+	  v[i][2] += dtfm*f[i][2];
+	}
       }
     }
   }
@@ -595,17 +640,29 @@ void FixNPH::final_integrate_respa(int ilevel)
 
     double **v = atom->v;
     double **f = atom->f;
+    double *rmass = atom->rmass;
     double *mass = atom->mass;
     int *type = atom->type;
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[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];
+	}
+      }
+    } else {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] += dtfm*f[i][0];
+	  v[i][1] += dtfm*f[i][1];
+	  v[i][2] += dtfm*f[i][2];
+	}
       }
     }
   }

src/fix_npt.cpp

@@ -262,8 +262,6 @@ int FixNPT::setmask()
 void FixNPT::init()
 {
   if (domain->triclinic) error->all("Cannot use fix npt with triclinic box");
-  if (atom->mass == NULL)
-    error->all("Cannot use fix npt without per-type mass defined");
 
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"deform") == 0) {
@@ -359,6 +357,7 @@ void FixNPT::setup(int vflag)
 void FixNPT::initial_integrate(int vflag)
 {
   int i;
+  double dtfm;
 
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
@@ -389,37 +388,59 @@ void FixNPT::initial_integrate(int vflag)
     dilation[i] = exp(dthalf*omega_dot[i]);
   }
 
-  // v update only for atoms in group
-
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  double dtfm;
+  if (rmass) {
-
+    if (which == NOBIAS) {
-  if (which == NOBIAS) {
+      for (i = 0; i < nlocal; i++) {
-    for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
-      if (mask[i] & groupbit) {
+	  dtfm = dtf / rmass[i];
-	dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
-	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
-	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
+      }
+    } else if (which == BIAS) {
+      for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / rmass[i];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
-  } else if (which == BIAS) {
+
-    for (i = 0; i < nlocal; i++) {
+  } else {
-      if (mask[i] & groupbit) {
+    if (which == NOBIAS) {
-	temperature->remove_bias(i,v[i]);
+      for (i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  dtfm = dtf / mass[type[i]];
-	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
-	temperature->restore_bias(i,v[i]);
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
+      }
+    } else if (which == BIAS) {
+      for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
   }
@@ -452,37 +473,60 @@ void FixNPT::initial_integrate(int vflag)
 void FixNPT::final_integrate()
 {
   int i;
-
+  double dtfm;
-  // v update only for atoms in group
 
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  double dtfm;
+  if (rmass) {
-
+    if (which == NOBIAS) {
-  if (which == NOBIAS) {
+      for (i = 0; i < nlocal; i++) {
-    for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
-      if (mask[i] & groupbit) {
+	  dtfm = dtf / rmass[i];
-	dtfm = dtf / mass[type[i]];
+	  v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
-	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	  v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
-	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	  v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
-	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	}
+      }
+    } else if (which == BIAS) {
+      for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / rmass[i];
+	  v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	  v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	  v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
-  } else if (which == BIAS) {
+
-    for (i = 0; i < nlocal; i++) {
+  } else {
-      if (mask[i] & groupbit) {
+    if (which == NOBIAS) {
-	temperature->remove_bias(i,v[i]);
+      for (i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	  dtfm = dtf / mass[type[i]];
-	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	  v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
-	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	  v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
-	temperature->restore_bias(i,v[i]);
+	  v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	}
+      }
+    } else if (which == BIAS) {
+      for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	  v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	  v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
   }
@@ -549,6 +593,7 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -592,12 +637,23 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
 
     // v update only for atoms in group
 
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  dtfm = dtf / rmass[i];
-	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
+      }
+    } else {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
       }
     }
 
@@ -609,24 +665,49 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
 
     // v update only for atoms in group
 
-    if (which == NOBIAS) {
+    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
+      if (which == NOBIAS) {
-	if (mask[i] & groupbit) {
+	for (int i = 0; i < nlocal; i++) {
-	  dtfm = dtf / mass[type[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];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[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];
+	    temperature->restore_bias(i,v[i]);
+	  }
 	}
       }
-    } else if (which == BIAS) {
+
-      for (int i = 0; i < nlocal; i++) {
+    } else {
-	if (mask[i] & groupbit) {
+      if (which == NOBIAS) {
-	  temperature->remove_bias(i,v[i]);
+	for (int i = 0; i < nlocal; i++) {
-	  dtfm = dtf / mass[type[i]];
+	  if (mask[i] & groupbit) {
-	  v[i][0] += dtfm*f[i][0];
+	    dtfm = dtf / mass[type[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];
-	  temperature->restore_bias(i,v[i]);
+	    v[i][2] += dtfm*f[i][2];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[type[i]];
+	    v[i][0] += dtfm*f[i][0];
+	    v[i][1] += dtfm*f[i][1];
+	    v[i][2] += dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	  }
 	}
       }
     }
@@ -668,30 +749,56 @@ void FixNPT::final_integrate_respa(int ilevel)
 
     double **v = atom->v;
     double **f = atom->f;
+    double *rmass = atom->rmass;
     double *mass = atom->mass;
     int *type = atom->type;
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
     if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-    if (which == NOBIAS) {
+    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
+      if (which == NOBIAS) {
-	if (mask[i] & groupbit) {
+	for (int i = 0; i < nlocal; i++) {
-	  dtfm = dtf / mass[type[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];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[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];
+	    temperature->restore_bias(i,v[i]);
+	  }
 	}
       }
-    } else if (which == BIAS) {
+
-      for (int i = 0; i < nlocal; i++) {
+    } else {
-	if (mask[i] & groupbit) {
+      if (which == NOBIAS) {
-	  temperature->remove_bias(i,v[i]);
+	for (int i = 0; i < nlocal; i++) {
-	  dtfm = dtf / mass[type[i]];
+	  if (mask[i] & groupbit) {
-	  v[i][0] += dtfm*f[i][0];
+	    dtfm = dtf / mass[type[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];
-	  temperature->restore_bias(i,v[i]);
+	    v[i][2] += dtfm*f[i][2];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[type[i]];
+	    v[i][0] += dtfm*f[i][0];
+	    v[i][1] += dtfm*f[i][1];
+	    v[i][2] += dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	  }
 	}
       }
     }

src/fix_npt_sphere.cpp

@@ -35,41 +35,61 @@ enum{NOBIAS,BIAS};
 FixNPTSphere::FixNPTSphere(LAMMPS *lmp, int narg, char **arg) :
   FixNPT(lmp, narg, arg)
 {
+  // error checks
+
   if (!atom->omega_flag || !atom->torque_flag)
     error->all("Fix npt/sphere requires atom attributes omega, torque");
-
+  if (!atom->radius_flag && !atom->avec->shape_type)
-  dttype = new double[atom->ntypes+1];
+    error->all("Fix npt/sphere requires atom attribute radius or shape");
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixNPTSphere::~FixNPTSphere()
-{
-  delete [] dttype;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNPTSphere::init()
 {
+  int i,itype;
+
+  // check that all particles are finite-size and spherical
+  // no point particles allowed
+
+  if (atom->radius_flag) {
+    double *radius = atom->radius;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+    if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	if (radius[i] == 0.0)
+	  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;
+    if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+    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");
+      }
+  }
+
   FixNPT::init();
-
-  if (!atom->shape)
-    error->all("Fix npt/sphere requires atom attribute shape");
-
-  double **shape = atom->shape;
-  for (int i = 1; i <= atom->ntypes; i++)
-    if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
-      error->all("Fix npt/sphere requires spherical particle shapes");
 }
 
-/* ----------------------------------------------------------------------
+/* ---------------------------------------------------------------------- */
-   1st half of Verlet update 
-------------------------------------------------------------------------- */
 
 void FixNPTSphere::initial_integrate(int vflag)
 {
-  int i;
+  int i,itype;
   double dtfm,dtirotate;
 
   double delta = update->ntimestep - update->beginstep;
@@ -102,37 +122,63 @@ void FixNPTSphere::initial_integrate(int vflag)
   }
   factor_rotate = exp(-dthalf*eta_dot);
 
-  // v update only for atoms in group
-
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
   double **omega = atom->omega;
   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;
 
-  if (which == NOBIAS) {
+  if (rmass) {
-    for (i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  dtfm = dtf / rmass[i];
-	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
+      }
+    } else if (which == BIAS) {
+      for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / rmass[i];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
-  } else if (which == BIAS) {
+
-    for (i = 0; i < nlocal; i++) {
+  } else {
-      if (mask[i] & groupbit) {
+    if (which == NOBIAS) {
-	temperature->remove_bias(i,v[i]);
+      for (i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  dtfm = dtf / mass[type[i]];
-	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
-	temperature->restore_bias(i,v[i]);
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	}
+      }
+    } else if (which == BIAS) {
+      for (i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
   }
@@ -151,24 +197,57 @@ void FixNPTSphere::initial_integrate(int vflag)
     }
   }
 
-  // recompute timesteps since dt may have changed or come via rRESPA
+  // set timestep here since dt may have changed or come via rRESPA
 
   double dtfrotate = dtf / INERTIA;
-  int ntypes = atom->ntypes;
-  double **shape = atom->shape;
-  for (int i = 1; i <= ntypes; i++)
-    dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
 
-  // update angular momentum by 1/2 step
+  // update omega for all particles
-  // update quaternion a full step via Richardson iteration
+  // d_omega/dt = torque / inertia
-  // returns new normalized quaternion
+  // 4 cases depending on radius vs shape and rmass vs mass
 
-  for (i = 0; i < nlocal; i++) {    
+  if (radius) {
-    if (mask[i] & groupbit) {
+    if (rmass) {
-      dtirotate = dttype[type[i]];
+      for (i = 0; i < nlocal; i++) {    
-      omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0];
+	if (mask[i] & groupbit) {
-      omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1];
+	  dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
-      omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2];
+	  omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0];
+	  omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1];
+	  omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2];
+	}
+      }
+    } else {
+      for (i = 0; i < nlocal; i++) {    
+	if (mask[i] & groupbit) {
+	  dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]);
+	  omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0];
+	  omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1];
+	  omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2];
+	}
+      }
+    }
+
+  } else {
+    if (rmass) {
+      for (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] = omega[i][0]*factor_rotate + dtirotate*torque[i][0];
+	  omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1];
+	  omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2];
+	}
+      }
+    } else {
+      for (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] = omega[i][0]*factor_rotate + dtirotate*torque[i][0];
+	  omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1];
+	  omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2];
+	}
+      }
     }
   }
 
@@ -179,68 +258,199 @@ void FixNPTSphere::initial_integrate(int vflag)
   if (kspace_flag) force->kspace->setup();
 }
 
-/* ----------------------------------------------------------------------
+/* ---------------------------------------------------------------------- */
-   2nd half of Verlet update 
-------------------------------------------------------------------------- */
 
 void FixNPTSphere::final_integrate()
 {
   int i,itype;
   double dtfm,dtirotate;
 
-  // v update only for atoms in group
-
   double **v = atom->v;
   double **f = atom->f;
   double **omega = atom->omega;
   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;
 
-  // recompute timesteps since dt may have changed or come via rRESPA
+  // set timestep here since dt may have changed or come via rRESPA
 
   double dtfrotate = dtf / INERTIA;
-  int ntypes = atom->ntypes;
-  double **shape = atom->shape;
-  for (int i = 1; i <= ntypes; i++)
-    dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
 
-  if (which == NOBIAS) {
+  // update v,omega for all particles
-    for (i = 0; i < nlocal; i++) {
+  // d_omega/dt = torque / inertia
-      if (mask[i] & groupbit) {
+  // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias
-	itype = type[i];
+
-	
+  if (radius) {
-	dtfm = dtf / mass[itype];
+    if (rmass) {
-	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+      if (which == NOBIAS) {
-	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	for (i = 0; i < nlocal; i++) {
-	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	  if (mask[i] & groupbit) {
-	
+	    dtfm = dtf / rmass[i];
-	dtirotate = dttype[itype];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
-	omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor_rotate;
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
-	omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor_rotate;
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
-	omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor_rotate;
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
+      } else if (which == BIAS) {
+	for (i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / rmass[i];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
+      }
+
+    } else {
+      if (which == NOBIAS) {
+	for (i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
+      } else if (which == BIAS) {
+	for (i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
       }
     }
 
-  } else if (which == BIAS) {
+  } else {
-    for (i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      if (which == NOBIAS) {
-	itype = type[i];
+	for (i = 0; i < nlocal; i++) {
-	
+	  if (mask[i] & groupbit) {
-	temperature->remove_bias(i,v[i]);
+	    itype = type[i];
-	dtfm = dtf / mass[itype];
+	    dtfm = dtf / rmass[i];
-	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
-	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
-	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
-	temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
+      } else if (which == BIAS) {
+	for (i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / rmass[i];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
+      }
 
-	dtirotate = dttype[itype];
+    } else {
-	omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor_rotate;
+      if (which == NOBIAS) {
-	omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor_rotate;
+	for (i = 0; i < nlocal; i++) {
-	omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor_rotate;
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    
+	    dtirotate = dtfrotate / 
+	      (shape[itype][0]*shape[itype][0]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
+      } else if (which == BIAS) {
+	for (i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	    v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	    v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / 
+	      (shape[itype][0]*shape[itype][0]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) *
+	      factor_rotate;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) *
+	      factor_rotate;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) *
+	      factor_rotate;
+	  }
+	}
       }
     }
   }

src/fix_npt_sphere.h

@@ -21,13 +21,12 @@ namespace LAMMPS_NS {
 class FixNPTSphere : public FixNPT {
  public:
   FixNPTSphere(class LAMMPS *, int, char **);
-  ~FixNPTSphere();
+  ~FixNPTSphere() {}
   void init();
   void initial_integrate(int);
   void final_integrate();
 
  private:
-  double *dttype;
   double factor_rotate;
 };
 

src/fix_nve_sphere.cpp

@@ -51,21 +51,14 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) :
     } else error->all("Illegal fix nve/sphere command");
   }
 
-  // error check
+  // error checks
 
   if (!atom->omega_flag || !atom->torque_flag)
     error->all("Fix nve/sphere requires atom attributes omega, torque");
+  if (!atom->radius_flag && !atom->avec->shape_type)
+    error->all("Fix nve/sphere requires atom attribute radius or shape");
   if (extra == DIPOLE && !atom->mu_flag)
     error->all("Fix nve/sphere requires atom attribute mu");
-
-  dttype = new double[atom->ntypes+1];
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixNVESphere::~FixNVESphere()
-{
-  delete [] dttype;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -84,23 +77,42 @@ int FixNVESphere::setmask()
 
 void FixNVESphere::init()
 {
-  dtv = update->dt;
+  int i,itype;
-  dtf = 0.5 * update->dt * force->ftm2v;
 
-  if (strcmp(update->integrate_style,"respa") == 0)
+  // check that all particles are finite-size and spherical
-    step_respa = ((Respa *) update->integrate)->step;
+  // no point particles allowed
 
-  if (atom->mass && !atom->shape)
+  if (atom->radius_flag) {
-    error->all("Fix nve/sphere requires atom attribute shape");
+    double *radius = atom->radius;
-  if (atom->rmass && !atom->radius_flag)
+    int *mask = atom->mask;
-    error->all("Fix nve/sphere requires atom attribute radius");
+    int nlocal = atom->nlocal;
+    if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  if (atom->mass) {
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	if (radius[i] == 0.0)
+	  error->one("Fix nve/sphere requires extended particles");
+      }
+
+  } else {
     double **shape = atom->shape;
-    for (int i = 1; i <= atom->ntypes; i++)
+    int *type = atom->type;
-      if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
+    int *mask = atom->mask;
-	error->all("Fix nve/sphere requires spherical particle shapes");
+    int nlocal = atom->nlocal;
+    if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+    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();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -119,58 +131,97 @@ void FixNVESphere::initial_integrate(int vflag)
   double *radius = atom->radius;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
-  int *mask = atom->mask;
+  double **shape = atom->shape;
   int *type = atom->type;
+  int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  // recompute timesteps since dt may have changed or come via rRESPA
+  // set timestep here since dt may have changed or come via rRESPA
 
-  dtfrotate = dtf / INERTIA;
+  double dtfrotate = dtf / INERTIA;
-  if (mass) {
-    double **shape = atom->shape;
-    int ntypes = atom->ntypes;
-    for (int i = 1; i <= ntypes; i++)
-      dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
-  }
 
   // update v,x,omega for all particles
   // d_omega/dt = torque / inertia
+  // 4 cases depending on radius vs shape and rmass vs mass
 
-  if (rmass) {
+  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 {
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	itype = type[i];
+	if (mask[i] & groupbit) {
-	dtfm = dtf / mass[itype];
+	  itype = type[i];
-	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 = dttype[itype];
+	  
-	omega[i][0] += dtirotate * torque[i][0];
+	  dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*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 / 
+	    (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];
+	}
       }
     }
   }
@@ -213,50 +264,85 @@ void FixNVESphere::final_integrate()
   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;
 
-  // recompute timesteps since dt may have changed or come via rRESPA
+  // set timestep here since dt may have changed or come via rRESPA
 
-  dtfrotate = dtf / INERTIA;
+  double dtfrotate = dtf / INERTIA;
-  if (mass) {
-    double **shape = atom->shape;
-    int ntypes = atom->ntypes;
-    for (int i = 1; i <= ntypes; i++)
-      dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
-  }
 
-  if (rmass) {
+  // update v,omega for all particles
-    for (int i = 0; i < nlocal; i++) {
+  // d_omega/dt = torque / inertia
-      if (mask[i] & groupbit) {
+  // 4 cases depending on radius vs shape and rmass vs mass
-	dtfm = dtf / rmass[i];
+
-	v[i][0] += dtfm * f[i][0];
+  if (radius) {
-	v[i][1] += dtfm * f[i][1];
+    if (rmass) {
-	v[i][2] += dtfm * f[i][2];
+      for (int i = 0; i < nlocal; i++) {
-	
+	if (mask[i] & groupbit) {
-	dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+	  dtfm = dtf / rmass[i];
-	omega[i][0] += dtirotate * torque[i][0];
+	  v[i][0] += dtfm * f[i][0];
-	omega[i][1] += dtirotate * torque[i][1];
+	  v[i][1] += dtfm * f[i][1];
-	omega[i][2] += dtirotate * torque[i][2];
+	  v[i][2] += dtfm * f[i][2];
+	  
+	  dtirotate = dtfrotate / (radius[i]*radius[i]*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 / (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 {
-    dtfrotate = dtf / INERTIA;
+    if (rmass) {
-    for (int i = 0; i < nlocal; i++) {
+      for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
+	if (mask[i] & groupbit) {
-	itype = type[i];
+	  itype = type[i];
-	dtfm = dtf / mass[itype];
+	  dtfm = dtf / rmass[i];
-	v[i][0] += dtfm * f[i][0];
+	  v[i][0] += dtfm * f[i][0];
-	v[i][1] += dtfm * f[i][1];
+	  v[i][1] += dtfm * f[i][1];
-	v[i][2] += dtfm * f[i][2];
+	  v[i][2] += dtfm * f[i][2];
-	
+	  
-	dtirotate = dttype[itype];
+	  dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
-	omega[i][0] += dtirotate * torque[i][0];
+	  omega[i][0] += dtirotate * torque[i][0];
-	omega[i][1] += dtirotate * torque[i][1];
+	  omega[i][1] += dtirotate * torque[i][1];
-	omega[i][2] += dtirotate * torque[i][2];
+	  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_nve_sphere.h

@@ -21,7 +21,7 @@ namespace LAMMPS_NS {
 class FixNVESphere : public FixNVE {
  public:
   FixNVESphere(class LAMMPS *, int, char **);
-  ~FixNVESphere();
+  ~FixNVESphere() {}
   int setmask();
   void init();
   void initial_integrate(int);
@@ -29,8 +29,6 @@ class FixNVESphere : public FixNVE {
 
  private:
   int extra;
-  double dtfrotate;
-  double *dttype;
 };
 
 }

src/fix_nvt.cpp

@@ -117,9 +117,6 @@ int FixNVT::setmask()
 
 void FixNVT::init()
 {
-  if (atom->mass == NULL)
-    error->all("Cannot use fix nvt without per-type mass defined");
-
   int icompute = modify->find_compute(id_temp);
   if (icompute < 0) error->all("Temp ID for fix nvt does not exist");
   temperature = modify->compute[icompute];
@@ -132,7 +129,6 @@ void FixNVT::init()
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
   dthalf = 0.5 * update->dt;
-
   drag_factor = 1.0 - (update->dt * t_freq * drag);
 
   if (strcmp(update->integrate_style,"respa") == 0) {
@@ -172,37 +168,70 @@ void FixNVT::initial_integrate(int vflag)
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  if (which == NOBIAS) {
+  if (rmass) {
-    for (int i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  dtfm = dtf / rmass[i];
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	x[i][0] += dtv * v[i][0];
+	  v[i][2] = v[i][2]*factor + 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];
+	}
+      }
+
+    } else if (which == BIAS) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / rmass[i];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	  x[i][0] += dtv * v[i][0];
+	  x[i][1] += dtv * v[i][1];
+	  x[i][2] += dtv * v[i][2];
+	}
       }
     }
 
-  } else if (which == BIAS) {
+  } else {
-    for (int i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	temperature->remove_bias(i,v[i]);
+	if (mask[i] & groupbit) {
-	dtfm = dtf / mass[type[i]];
+	  dtfm = dtf / mass[type[i]];
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(i,v[i]);
+	  x[i][0] += dtv * v[i][0];
-	x[i][0] += dtv * v[i][0];
+	  x[i][1] += dtv * v[i][1];
-	x[i][1] += dtv * v[i][1];
+	  x[i][2] += dtv * v[i][2];
-	x[i][2] += dtv * v[i][2];
+	}
+      }
+
+    } else if (which == BIAS) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	  x[i][0] += dtv * v[i][0];
+	  x[i][1] += dtv * v[i][1];
+	  x[i][2] += dtv * v[i][2];
+	}
       }
     }
   }
@@ -218,31 +247,58 @@ void FixNVT::final_integrate()
 
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  if (which == NOBIAS) {
+  if (rmass) {
-    for (int i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]] * factor;
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  dtfm = dtf / rmass[i] * factor;
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	}
+      }
+      
+    } else if (which == BIAS) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / rmass[i] * factor;
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
 
-  } else if (which == BIAS) {
+  } else {
-    for (int i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	temperature->remove_bias(i,v[i]);
+	if (mask[i] & groupbit) {
-	dtfm = dtf / mass[type[i]] * factor;
+	  dtfm = dtf / mass[type[i]] * factor;
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(i,v[i]);
+	}
+      }
+      
+    } else if (which == BIAS) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / mass[type[i]] * factor;
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
   }
@@ -276,6 +332,7 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -300,25 +357,51 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
     factor = exp(-dthalf*eta_dot);
   } else factor = 1.0;
 
-  if (which == NOBIAS) {
+  if (rmass) {
-    for (int i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[i]];
+	if (mask[i] & groupbit) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  dtfm = dtf / rmass[i];
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	}
+      }
+
+    } else if (which == BIAS) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / rmass[i];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
 
-  } else if (which == BIAS) {
+  } else {
-    for (int i = 0; i < nlocal; i++) {
+    if (which == NOBIAS) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	temperature->remove_bias(i,v[i]);
+	if (mask[i] & groupbit) {
-	dtfm = dtf / mass[type[i]];
+	  dtfm = dtf / mass[type[i]];
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(i,v[i]);
+	}
+      }
+
+    } else if (which == BIAS) {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  temperature->remove_bias(i,v[i]);
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	  v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	  temperature->restore_bias(i,v[i]);
+	}
       }
     }
   }
@@ -351,18 +434,31 @@ void FixNVT::final_integrate_respa(int ilevel)
   else {
     double **v = atom->v;
     double **f = atom->f;
+    double *rmass = atom->rmass;
     double *mass = atom->mass;
     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 (rmass) {
-      if (mask[i] & groupbit) {
+      for (int i = 0; i < nlocal; i++) {
-	dtfm = dtf / mass[type[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];
+	}
+      }
+      
+    } else {
+      for (int i = 0; i < nlocal; i++) {
+	if (mask[i] & groupbit) {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] += dtfm*f[i][0];
+	  v[i][1] += dtfm*f[i][1];
+	  v[i][2] += dtfm*f[i][2];
+	}
       }
     }
   }
@@ -440,7 +536,6 @@ void FixNVT::reset_dt()
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
   dthalf = 0.5 * update->dt;
-
   drag_factor = 1.0 - (update->dt * t_freq * drag);
 }
 

src/fix_nvt_sphere.cpp

@@ -36,32 +36,54 @@ enum{NOBIAS,BIAS};
 FixNVTSphere::FixNVTSphere(LAMMPS *lmp, int narg, char **arg) :
   FixNVT(lmp, narg, arg)
 {
+  // error checks
+
   if (!atom->omega_flag || !atom->torque_flag)
     error->all("Fix nvt/sphere requires atom attributes omega, torque");
-
+  if (!atom->radius_flag && !atom->avec->shape_type)
-  dttype = new double[atom->ntypes+1];
+    error->all("Fix nvt/sphere requires atom attribute radius or shape");
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixNVTSphere::~FixNVTSphere()
-{
-  delete [] dttype;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNVTSphere::init()
 {
+  int i,itype;
+
+  // check that all particles are finite-size and spherical
+  // no point particles allowed
+
+  if (atom->radius_flag) {
+    double *radius = atom->radius;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+    if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	if (radius[i] == 0.0)
+	  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;
+    if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+    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");
+      }
+  }
+
   FixNVT::init();
-
-  if (!atom->shape)
-    error->all("Fix nvt/sphere requires atom attribute shape");
-
-  double **shape = atom->shape;
-  for (int i = 1; i <= atom->ntypes; i++)
-    if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
-      error->all("Fix nvt/sphere requires spherical particle shapes");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -90,59 +112,188 @@ void FixNVTSphere::initial_integrate(int vflag)
   double **f = atom->f;
   double **omega = atom->omega;
   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;
 
-  // recompute timesteps since dt may have changed or come via rRESPA
+  // set timestep here since dt may have changed or come via rRESPA
 
   double dtfrotate = dtf / INERTIA;
-  int ntypes = atom->ntypes;
-  double **shape = atom->shape;
-  for (int i = 1; i <= ntypes; i++)
-    dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
 
-  if (which == NOBIAS) {
+  // update v,x,omega for all particles
-    for (int i = 0; i < nlocal; i++) {
+  // d_omega/dt = torque / inertia
-      if (mask[i] & groupbit) {
+  // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias
-	itype = type[i];
 
-	dtfm = dtf / mass[itype];
+  if (radius) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+    if (rmass) {
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+      if (which == NOBIAS) {
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	for (int i = 0; i < nlocal; i++) {
-	x[i][0] += dtv * v[i][0];
+	  if (mask[i] & groupbit) {
-	x[i][1] += dtv * v[i][1];
+	    dtfm = dtf / rmass[i];
-	x[i][2] += dtv * v[i][2];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	dtirotate = dttype[itype];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    x[i][0] += dtv * v[i][0];
-	omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    x[i][1] += dtv * v[i][1];
-	omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	    x[i][2] += dtv * v[i][2];
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+	    omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / rmass[i];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    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]*rmass[i]);
+	    omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
+      }
+
+    } else {
+      if (which == NOBIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + 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] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    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] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
       }
     }
 
-  } else if (which == BIAS) {
+  } else {
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      if (which == NOBIAS) {
-	itype = type[i];
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    dtfm = dtf / rmass[i];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + 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] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / rmass[i];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    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] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
+      }
 
-	temperature->remove_bias(i,v[i]);
+    } else {
-	dtfm = dtf / mass[itype];
+      if (which == NOBIAS) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	for (int i = 0; i < nlocal; i++) {
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  if (mask[i] & groupbit) {
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    itype = type[i];
-	temperature->restore_bias(i,v[i]);
+	    dtfm = dtf / mass[itype];
-	x[i][0] += dtv * v[i][0];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	x[i][1] += dtv * v[i][1];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	x[i][2] += dtv * v[i][2];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	
+	    x[i][0] += dtv * v[i][0];
-	dtirotate = dttype[itype];
+	    x[i][1] += dtv * v[i][1];
-	omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    x[i][2] += dtv * v[i][2];
-	omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    
-	omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	    dtirotate = dtfrotate / 
+	      (shape[itype][0]*shape[itype][0]*mass[itype]);
+	    omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
+      } else if (which == BIAS) {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[itype];
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    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] = omega[i][0]*factor + dtirotate*torque[i][0];
+	    omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1];
+	    omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2];
+	  }
+	}
       }
     }
   }
@@ -161,53 +312,164 @@ void FixNVTSphere::final_integrate()
   double **f = atom->f;
   double **omega = atom->omega;
   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;
 
-  // recompute timesteps since dt may have changed or come via rRESPA
+  // set timestep here since dt may have changed or come via rRESPA
 
   double dtfrotate = dtf / INERTIA;
-  int ntypes = atom->ntypes;
-  double **shape = atom->shape;
-  for (int i = 1; i <= ntypes; i++)
-    dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
 
-  if (which == NOBIAS) {
+  // update v,omega for all particles
-    for (int i = 0; i < nlocal; i++) {
+  // d_omega/dt = torque / inertia
-      if (mask[i] & groupbit) {
+  // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias
-	itype = type[i];
 
-	dtfm = dtf / mass[itype] * factor;
+  if (radius) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+    if (rmass) {
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+      if (which == NOBIAS) {
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    dtfm = dtf / rmass[i] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
+      } else {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / rmass[i] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
+      }
 
-	dtirotate = dttype[itype];
+    } else {
-	omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+      if (which == NOBIAS) {
-	omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	for (int i = 0; i < nlocal; i++) {
-	omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    dtfm = dtf / mass[itype] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
+      } else {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[itype] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
       }
     }
 
   } else {
-    for (int i = 0; i < nlocal; i++) {
+    if (rmass) {
-      if (mask[i] & groupbit) {
+      if (which == NOBIAS) {
-	itype = type[i];
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    dtfm = dtf / rmass[i] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    
+	    dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
+      } else {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / rmass[i] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
+      }
 
-	temperature->remove_bias(i,v[i]);
+    } else {
-	dtfm = dtf / mass[itype] * factor;
+      if (which == NOBIAS) {
-	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	for (int i = 0; i < nlocal; i++) {
-	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	  if (mask[i] & groupbit) {
-	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    itype = type[i];
-	temperature->restore_bias(i,v[i]);
+	    dtfm = dtf / mass[itype] * factor;
-	
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-	dtirotate = dttype[itype];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-	omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    
-	omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	    dtirotate = dtfrotate / 
+	      (shape[itype][0]*shape[itype][0]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
+      } else {
+	for (int i = 0; i < nlocal; i++) {
+	  if (mask[i] & groupbit) {
+	    itype = type[i];
+	    temperature->remove_bias(i,v[i]);
+	    dtfm = dtf / mass[itype] * factor;
+	    v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	    v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	    v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	    temperature->restore_bias(i,v[i]);
+	    
+	    dtirotate = dtfrotate / 
+	      (shape[itype][0]*shape[itype][0]*mass[itype]);
+	    omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor;
+	    omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor;
+	    omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor;
+	  }
+	}
       }
     }
   }

src/fix_nvt_sphere.h

@@ -21,13 +21,10 @@ namespace LAMMPS_NS {
 class FixNVTSphere : public FixNVT {
  public:
   FixNVTSphere(class LAMMPS *, int, char **);
-  ~FixNVTSphere();
+  ~FixNVTSphere() {}
   void init();
   void initial_integrate(int);
   void final_integrate();
-
- private:
-  double *dttype;
 };
 
 }

src/fix_press_berendsen.cpp

@@ -234,15 +234,14 @@ void FixPressBerendsen::init()
 {
   if (domain->triclinic) 
     error->all("Cannot use fix press/berendsen with triclinic box");
-  if (atom->mass == NULL)
-    error->all("Cannot use fix press/berendsen without per-type mass defined");
 
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"deform") == 0) {
       int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
       if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || 
 	  (p_flag[2] && dimflag[2]))
-	error->all("Cannot use fix press/berendsen and fix deform on same dimension");
+	error->all("Cannot use fix press/berendsen and "
+		   "fix deform on same dimension");
     }
 
   // set temperature and pressure ptrs

src/fix_rigid.cpp

@@ -368,9 +368,6 @@ void FixRigid::init()
 
   triclinic = domain->triclinic;
 
-  if (atom->mass == NULL)
-    error->all("Cannot use fix rigid without per-type mass defined");
-
   // warn if more than one rigid fix
 
   int count = 0;
@@ -401,10 +398,11 @@ void FixRigid::init()
 
   // compute masstotal & center-of-mass of each rigid body
   
-  int *type = atom->type;
-  int *image = atom->image;
-  double *mass = atom->mass;
   double **x = atom->x;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *image = atom->image;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
   
   double xprd = domain->xprd;
@@ -426,7 +424,8 @@ void FixRigid::init()
     xbox = (image[i] & 1023) - 512;
     ybox = (image[i] >> 10 & 1023) - 512;
     zbox = (image[i] >> 20) - 512;
-    massone = mass[type[i]];
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
     
     if (triclinic == 0) {
       xunwrap = x[i][0] + xbox*xprd;
@@ -486,7 +485,8 @@ void FixRigid::init()
     dx = xunwrap - xcm[ibody][0];
     dy = yunwrap - xcm[ibody][1];
     dz = zunwrap - xcm[ibody][2];
-    massone = mass[type[i]];
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
     
     sum[ibody][0] += massone * (dy*dy + dz*dz);
     sum[ibody][1] += massone * (dx*dx + dz*dz);
@@ -613,7 +613,8 @@ void FixRigid::init()
   for (i = 0; i < nlocal; i++) {
     if (body[i] < 0) continue;
     ibody = body[i];
-    massone = mass[type[i]];
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
 
     sum[ibody][0] += massone * 
       (displace[i][1]*displace[i][1] + displace[i][2]*displace[i][2]);
@@ -669,6 +670,7 @@ void FixRigid::setup(int vflag)
   int *type = atom->type;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass;
   double *mass = atom->mass;
   int nlocal = atom->nlocal;
 
@@ -679,7 +681,9 @@ void FixRigid::setup(int vflag)
   for (i = 0; i < nlocal; i++) {
     if (body[i] < 0) continue;
     ibody = body[i];
-    massone = mass[type[i]];
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
+
     sum[ibody][0] += v[i][0] * massone;
     sum[ibody][1] += v[i][1] * massone;
     sum[ibody][2] += v[i][2] * massone;
@@ -739,7 +743,9 @@ void FixRigid::setup(int vflag)
     dy = yunwrap - xcm[ibody][1];
     dz = zunwrap - xcm[ibody][2];
     
-    massone = mass[type[i]];
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
+
     sum[ibody][0] += dy * massone*v[i][2] - dz * massone*v[i][1];
     sum[ibody][1] += dz * massone*v[i][0] - dx * massone*v[i][2];
     sum[ibody][2] += dx * massone*v[i][1] - dy * massone*v[i][0];
@@ -1413,6 +1419,7 @@ void FixRigid::set_xv()
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
+  double *rmass = atom->rmass; 
   double *mass = atom->mass; 
   int *type = atom->type;
   int nlocal = atom->nlocal;
@@ -1495,7 +1502,9 @@ void FixRigid::set_xv()
     // assume per-atom contribution is due to constraint force on that atom
 
     if (evflag) {
-      massone = mass[type[i]];
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+
       fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
       fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
       fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; 
@@ -1526,10 +1535,11 @@ void FixRigid::set_v()
   double xy,xz,yz;
   double vr[6];
 
-  double *mass = atom->mass; 
   double **f = atom->f;
   double **v = atom->v;
   double **x = atom->x;
+  double *rmass = atom->rmass; 
+  double *mass = atom->mass; 
   int *type = atom->type;
   int *image = atom->image;
   int nlocal = atom->nlocal;
@@ -1578,7 +1588,9 @@ void FixRigid::set_v()
     // assume per-atom contribution is due to constraint force on that atom
 
     if (evflag) {
-      massone = mass[type[i]];
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+
       fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
       fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
       fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; 

src/fix_temp_berendsen.cpp

@@ -87,9 +87,6 @@ int FixTempBerendsen::setmask()
 
 void FixTempBerendsen::init()
 {
-  if (atom->mass == NULL)
-    error->all("Cannot use fix temp/berendsen without per-type mass defined");
-
   int icompute = modify->find_compute(id_temp);
   if (icompute < 0)
     error->all("Temp ID for fix temp/berendsen does not exist");
@@ -127,7 +124,6 @@ void FixTempBerendsen::end_of_step()
 	v[i][2] *= lamda;
       }
     }
-
   } else if (which == BIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {

src/variable.cpp

@@ -1424,13 +1424,18 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
 			     double *argstack, int &nargstack)
 {
   // parse contents for arg1,arg2,arg3 separated by commas
-
+  // ptr1,ptr2 = location of 1st and 2nd comma, NULL if none
-  char *ptr1 = strchr(contents,',');
-  if (ptr1) *ptr1 = '\0';
-  char *ptr2 = strchr(ptr1+1,',');
-  if (ptr2) *ptr2 = '\0';
 
   char *arg1,*arg2,*arg3;
+  char *ptr1,*ptr2;
+
+  ptr1 = strchr(contents,',');
+  if (ptr1) {
+    *ptr1 = '\0';
+    ptr2 = strchr(ptr1+1,',');
+    if (ptr2) *ptr2 = '\0';
+  } else ptr2 = NULL;
+
   int n = strlen(contents) + 1;
   arg1 = new char[n];
   strcpy(arg1,contents);