diff -Naur lammps-9Apr08/src/ASPHERE/fix_npt_asphere.cpp lammps-10Apr08/src/ASPHERE/fix_npt_asphere.cpp
--- lammps-9Apr08/src/ASPHERE/fix_npt_asphere.cpp	2008-03-20 09:32:33.000000000 -0600
+++ lammps-10Apr08/src/ASPHERE/fix_npt_asphere.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -44,14 +44,6 @@
 	       "quat, angmom, torque, shape");
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixNPTAsphere::init()
-{
-  FixNPT::init();
-  dtq = 0.5 * update->dt;
-}
-
 /* ----------------------------------------------------------------------
    1st half of Verlet update 
 ------------------------------------------------------------------------- */
@@ -61,6 +53,8 @@
   int i;
   double dtfm;
 
+  dtq = 0.5 * dtv;
+
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
 
@@ -249,14 +243,6 @@
   }
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixNPTAsphere::reset_dt()
-{
-  FixNPT::reset_dt();
-  dtq = 0.5 * update->dt;
-}
-
 /* ----------------------------------------------------------------------
    Richardson iteration to update quaternion accurately
 ------------------------------------------------------------------------- */
diff -Naur lammps-9Apr08/src/ASPHERE/fix_npt_asphere.h lammps-10Apr08/src/ASPHERE/fix_npt_asphere.h
--- lammps-9Apr08/src/ASPHERE/fix_npt_asphere.h	2008-03-20 08:53:21.000000000 -0600
+++ lammps-10Apr08/src/ASPHERE/fix_npt_asphere.h	2008-04-08 16:55:40.000000000 -0600
@@ -21,10 +21,8 @@
 class FixNPTAsphere : public FixNPT {
  public:
   FixNPTAsphere(class LAMMPS *, int, char **);
-  void init();
   void initial_integrate(int);
   void final_integrate();
-  void reset_dt();
 
  private:
   double dtq;
diff -Naur lammps-9Apr08/src/ASPHERE/fix_nve_asphere.cpp lammps-10Apr08/src/ASPHERE/fix_nve_asphere.cpp
--- lammps-9Apr08/src/ASPHERE/fix_nve_asphere.cpp	2008-03-20 08:53:21.000000000 -0600
+++ lammps-10Apr08/src/ASPHERE/fix_nve_asphere.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -66,6 +66,8 @@
 {
   double dtfm;
 
+  dtq = 0.5 * dtv;
+
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
diff -Naur lammps-9Apr08/src/ASPHERE/fix_nve_asphere.h lammps-10Apr08/src/ASPHERE/fix_nve_asphere.h
--- lammps-9Apr08/src/ASPHERE/fix_nve_asphere.h	2008-03-20 08:53:21.000000000 -0600
+++ lammps-10Apr08/src/ASPHERE/fix_nve_asphere.h	2008-04-08 16:55:40.000000000 -0600
@@ -27,6 +27,7 @@
   void final_integrate();
 
  private:
+  double dtq;
   double **inertia;
 
   void richardson(double *, double *, double *);
diff -Naur lammps-9Apr08/src/ASPHERE/fix_nvt_asphere.cpp lammps-10Apr08/src/ASPHERE/fix_nvt_asphere.cpp
--- lammps-9Apr08/src/ASPHERE/fix_nvt_asphere.cpp	2008-03-20 09:32:33.000000000 -0600
+++ lammps-10Apr08/src/ASPHERE/fix_nvt_asphere.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -47,18 +47,12 @@
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVTAsphere::init()
-{
-  FixNVT::init();
-  dtq = 0.5 * update->dt;
-}
-
-/* ---------------------------------------------------------------------- */
-
 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);
@@ -197,14 +191,6 @@
   eta_dot *= drag_factor;
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixNVTAsphere::reset_dt()
-{
-  FixNVT::reset_dt();
-  dtq = 0.5 * update->dt;
-}
-
 /* ----------------------------------------------------------------------
    Richardson iteration to update quaternion accurately
 ------------------------------------------------------------------------- */
diff -Naur lammps-9Apr08/src/ASPHERE/fix_nvt_asphere.h lammps-10Apr08/src/ASPHERE/fix_nvt_asphere.h
--- lammps-9Apr08/src/ASPHERE/fix_nvt_asphere.h	2008-03-20 08:53:21.000000000 -0600
+++ lammps-10Apr08/src/ASPHERE/fix_nvt_asphere.h	2008-04-08 16:55:40.000000000 -0600
@@ -21,10 +21,8 @@
 class FixNVTAsphere : public FixNVT {
  public:
   FixNVTAsphere(class LAMMPS *, int, char **);
-  void init();
   void initial_integrate(int);
   void final_integrate();
-  void reset_dt();
 
  private:
   double dtq;
diff -Naur lammps-9Apr08/src/atom.cpp lammps-10Apr08/src/atom.cpp
--- lammps-9Apr08/src/atom.cpp	2008-04-08 07:58:15.000000000 -0600
+++ lammps-10Apr08/src/atom.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -601,26 +601,29 @@
 }
 
 /* ----------------------------------------------------------------------
-   check if atom tags are consecutive from 1 to Natoms
-   return 0 if any tag <= 0 or maxtag != Natoms
-   return 1 if OK (doesn't actually check if all tag values are used)
+   check that atom IDs span range from 1 to Natoms
+   return 0 if mintag != 1 or maxtag != Natoms
+   return 1 if OK
+   doesn't actually check if all tag values are used
 ------------------------------------------------------------------------- */
 
 int Atom::tag_consecutive()
 {
-  // check[0] = flagged if any tag <= 0
-  // check[1] = max tag of any atom
-
   int check[2],check_all[2];
-
   check[0] = check[1] = 0;
+
+  int idmin = static_cast<int> (natoms);
+  int idmax = 0;
+  
   for (int i = 0; i < nlocal; i++) {
-    if (tag[i] <= 0) check[0] = 1;
-    if (tag[i] > check[1]) check[1] = tag[i];
+    idmin = MIN(idmin,tag[i]);
+    idmax = MAX(idmax,tag[i]);
   }
-  MPI_Allreduce(check,check_all,2,MPI_INT,MPI_MAX,world);
+  int idminall,idmaxall;
+  MPI_Allreduce(&idmin,&idminall,1,MPI_INT,MPI_MIN,world);
+  MPI_Allreduce(&idmax,&idmaxall,1,MPI_INT,MPI_MAX,world);
 
-  if (check_all[0] || check_all[1] != static_cast<int> (natoms)) return 0;
+  if (idminall != 1 || idmaxall != static_cast<int> (natoms)) return 0;
   return 1;
 }
 
diff -Naur lammps-9Apr08/src/fix_add_force.cpp lammps-10Apr08/src/fix_add_force.cpp
--- lammps-9Apr08/src/fix_add_force.cpp	2008-03-17 14:49:18.000000000 -0600
+++ lammps-10Apr08/src/fix_add_force.cpp	2008-04-09 10:41:47.000000000 -0600
@@ -27,6 +27,12 @@
   Fix(lmp, narg, arg)
 {
   if (narg != 6) error->all("Illegal fix addforce command");
+
+  vector_flag = 1;
+  size_vector = 3;
+  scalar_vector_freq = 1;
+  extvector = 1;
+
   xvalue = atof(arg[3]);
   yvalue = atof(arg[4]);
   zvalue = atof(arg[5]);
diff -Naur lammps-9Apr08/src/fix_ave_force.cpp lammps-10Apr08/src/fix_ave_force.cpp
--- lammps-9Apr08/src/fix_ave_force.cpp	2008-03-17 14:49:18.000000000 -0600
+++ lammps-10Apr08/src/fix_ave_force.cpp	2008-04-09 10:41:47.000000000 -0600
@@ -29,6 +29,11 @@
 {
   if (narg != 6) error->all("Illegal fix aveforce command");
 
+  vector_flag = 1;
+  size_vector = 3;
+  scalar_vector_freq = 1;
+  extvector = 1;
+
   xflag = yflag = zflag = 1;
   if (strcmp(arg[3],"NULL") == 0) xflag = 0;
   else xvalue = atof(arg[3]);
diff -Naur lammps-9Apr08/src/fix_npt.h lammps-10Apr08/src/fix_npt.h
--- lammps-9Apr08/src/fix_npt.h	2008-03-19 15:43:26.000000000 -0600
+++ lammps-10Apr08/src/fix_npt.h	2008-04-08 16:55:40.000000000 -0600
@@ -25,7 +25,7 @@
   int setmask();
   virtual void init();
   void setup(int);
-  void initial_integrate(int);
+  virtual void initial_integrate(int);
   virtual void final_integrate();
   void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
@@ -33,7 +33,7 @@
   void write_restart(FILE *);
   void restart(char *);
   int modify_param(int, char **);
-  virtual void reset_dt();
+  void reset_dt();
 
  protected:
   int dimension,which;
diff -Naur lammps-9Apr08/src/fix_npt_sphere.cpp lammps-10Apr08/src/fix_npt_sphere.cpp
--- lammps-9Apr08/src/fix_npt_sphere.cpp	2008-03-20 19:04:20.000000000 -0600
+++ lammps-10Apr08/src/fix_npt_sphere.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -53,18 +53,14 @@
 void FixNPTSphere::init()
 {
   FixNPT::init();
-  dtfrotate = dtf / INERTIA;
 
   if (!atom->shape)
     error->all("Fix npt/sphere requires atom attribute shape");
 
-  double *mass = atom->mass;
   double **shape = atom->shape;
-  for (int i = 1; i <= atom->ntypes; i++) {
+  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");
-    dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]);
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -155,6 +151,14 @@
     }
   }
 
+  // recompute timesteps 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 / (0.25*shape[i][0]*shape[i][0]*mass[i]);
+
   // update angular momentum by 1/2 step
   // update quaternion a full step via Richardson iteration
   // returns new normalized quaternion
@@ -196,6 +200,14 @@
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // recompute timesteps 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 / (0.25*shape[i][0]*shape[i][0]*mass[i]);
+
   if (which == NOBIAS) {
     for (i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -268,11 +280,3 @@
     omega_dot[i] *= drag_factor;
   }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixNPTSphere::reset_dt()
-{
-  FixNPT::reset_dt();
-  dtfrotate = dtf / INERTIA;
-}
diff -Naur lammps-9Apr08/src/fix_npt_sphere.h lammps-10Apr08/src/fix_npt_sphere.h
--- lammps-9Apr08/src/fix_npt_sphere.h	2008-03-20 19:04:20.000000000 -0600
+++ lammps-10Apr08/src/fix_npt_sphere.h	2008-04-08 16:55:40.000000000 -0600
@@ -25,10 +25,8 @@
   void init();
   void initial_integrate(int);
   void final_integrate();
-  void reset_dt();
 
  private:
-  double dtfrotate;
   double *dttype;
   double factor_rotate;
 };
diff -Naur lammps-9Apr08/src/fix_nve.cpp lammps-10Apr08/src/fix_nve.cpp
--- lammps-9Apr08/src/fix_nve.cpp	2008-04-04 17:18:14.000000000 -0600
+++ lammps-10Apr08/src/fix_nve.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -50,7 +50,6 @@
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
 
   if (strcmp(update->integrate_style,"respa") == 0)
     step_respa = ((Respa *) update->integrate)->step;
@@ -169,5 +168,4 @@
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
 }
diff -Naur lammps-9Apr08/src/fix_nve.h lammps-10Apr08/src/fix_nve.h
--- lammps-9Apr08/src/fix_nve.h	2008-01-09 14:56:57.000000000 -0700
+++ lammps-10Apr08/src/fix_nve.h	2008-04-08 16:55:40.000000000 -0600
@@ -30,7 +30,7 @@
   void reset_dt();
 
  protected:
-  double dtv,dtf,dtq;
+  double dtv,dtf;
   double *step_respa;
   int mass_require;
 };
diff -Naur lammps-9Apr08/src/fix_nve_sphere.cpp lammps-10Apr08/src/fix_nve_sphere.cpp
--- lammps-9Apr08/src/fix_nve_sphere.cpp	2008-04-04 17:17:45.000000000 -0600
+++ lammps-10Apr08/src/fix_nve_sphere.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -31,7 +31,7 @@
 /* ---------------------------------------------------------------------- */
 
 FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
+  FixNVE(lmp, narg, arg)
 {
   if (narg < 3) error->all("Illegal fix nve/sphere command");
 
@@ -86,7 +86,6 @@
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtfrotate = dtf / INERTIA;
 
   if (strcmp(update->integrate_style,"respa") == 0)
     step_respa = ((Respa *) update->integrate)->step;
@@ -97,13 +96,10 @@
     error->all("Fix nve/sphere requires atom attributes radius, rmass");
 
   if (atom->mass) {
-    double *mass = atom->mass;
     double **shape = atom->shape;
-    for (int i = 1; i <= atom->ntypes; i++) {
+    for (int i = 1; i <= atom->ntypes; i++)
       if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
 	error->all("Fix nve/sphere requires spherical particle shapes");
-      dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]);
-    }
   }
 }
 
@@ -128,6 +124,16 @@
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // recompute timesteps since dt may have changed or come via rRESPA
+
+  dtfrotate = dtf / INERTIA;
+  if (mass) {
+    double **shape = atom->shape;
+    int ntypes = atom->ntypes;
+    for (int i = 1; i <= ntypes; i++)
+      dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]);
+  }
+
   // update v,x,omega for all particles
   // d_omega/dt = torque / inertia
 
@@ -212,6 +218,16 @@
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // recompute timesteps since dt may have changed or come via rRESPA
+
+  dtfrotate = dtf / INERTIA;
+  if (mass) {
+    double **shape = atom->shape;
+    int ntypes = atom->ntypes;
+    for (int i = 1; i <= ntypes; i++)
+      dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]);
+  }
+
   if (mass) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -229,6 +245,7 @@
     }
     
   } else {
+    dtfrotate = dtf / INERTIA;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
 	dtfm = dtf / rmass[i];
@@ -244,35 +261,3 @@
     }
   }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixNVESphere::initial_integrate_respa(int vflag, int ilevel, int flag)
-{
-  if (flag) return;             // only used by NPT,NPH
-
-  dtv = step_respa[ilevel];
-  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
-  dtfrotate = dtf / INERTIA;
-
-  if (ilevel == 0) initial_integrate(vflag);
-  else final_integrate();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixNVESphere::final_integrate_respa(int ilevel)
-{
-  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
-  dtfrotate = dtf / INERTIA;
-  final_integrate();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixNVESphere::reset_dt()
-{
-  dtv = update->dt;
-  dtf = 0.5 * update->dt * force->ftm2v;
-  dtfrotate = dtf / INERTIA;
-}
diff -Naur lammps-9Apr08/src/fix_nve_sphere.h lammps-10Apr08/src/fix_nve_sphere.h
--- lammps-9Apr08/src/fix_nve_sphere.h	2008-03-17 17:37:19.000000000 -0600
+++ lammps-10Apr08/src/fix_nve_sphere.h	2008-04-08 16:55:40.000000000 -0600
@@ -14,11 +14,11 @@
 #ifndef FIX_NVE_SPHERE_H
 #define FIX_NVE_SPHERE_H
 
-#include "fix.h"
+#include "fix_nve.h"
 
 namespace LAMMPS_NS {
 
-class FixNVESphere : public Fix {
+class FixNVESphere : public FixNVE {
  public:
   FixNVESphere(class LAMMPS *, int, char **);
   ~FixNVESphere();
@@ -26,14 +26,10 @@
   void init();
   void initial_integrate(int);
   void final_integrate();
-  void initial_integrate_respa(int, int, int);
-  void final_integrate_respa(int);
-  void reset_dt();
 
  private:
   int extra;
-  double dtv,dtf,dtfrotate;
-  double *step_respa;
+  double dtfrotate;
   double *dttype;
 };
 
diff -Naur lammps-9Apr08/src/fix_nvt.h lammps-10Apr08/src/fix_nvt.h
--- lammps-9Apr08/src/fix_nvt.h	2008-03-19 15:43:26.000000000 -0600
+++ lammps-10Apr08/src/fix_nvt.h	2008-04-08 16:55:40.000000000 -0600
@@ -34,7 +34,7 @@
   void restart(char *);
   int modify_param(int, char **);
   void reset_target(double);
-  virtual void reset_dt();
+  void reset_dt();
 
  protected:
   int which;
diff -Naur lammps-9Apr08/src/fix_nvt_sphere.cpp lammps-10Apr08/src/fix_nvt_sphere.cpp
--- lammps-9Apr08/src/fix_nvt_sphere.cpp	2008-03-20 19:04:20.000000000 -0600
+++ lammps-10Apr08/src/fix_nvt_sphere.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -54,18 +54,14 @@
 void FixNVTSphere::init()
 {
   FixNVT::init();
-  dtfrotate = dtf / INERTIA;
 
   if (!atom->shape)
     error->all("Fix nvt/sphere requires atom attribute shape");
 
-  double *mass = atom->mass;
   double **shape = atom->shape;
-  for (int i = 1; i <= atom->ntypes; i++) {
+  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");
-    dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]);
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -101,6 +97,14 @@
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // recompute timesteps 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 / (0.25*shape[i][0]*shape[i][0]*mass[i]);
+
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -164,6 +168,14 @@
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
+  // recompute timesteps 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 / (0.25*shape[i][0]*shape[i][0]*mass[i]);
+
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -211,11 +223,3 @@
   eta_dot += f_eta*dthalf;
   eta_dot *= drag_factor;
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixNVTSphere::reset_dt()
-{
-  FixNVT::reset_dt();
-  dtfrotate = dtf / INERTIA;
-}
diff -Naur lammps-9Apr08/src/fix_nvt_sphere.h lammps-10Apr08/src/fix_nvt_sphere.h
--- lammps-9Apr08/src/fix_nvt_sphere.h	2008-03-20 19:04:20.000000000 -0600
+++ lammps-10Apr08/src/fix_nvt_sphere.h	2008-04-08 16:55:40.000000000 -0600
@@ -25,10 +25,8 @@
   void init();
   void initial_integrate(int);
   void final_integrate();
-  void reset_dt();
 
  private:
-  double dtfrotate;
   double *dttype;
 };
 
diff -Naur lammps-9Apr08/src/velocity.cpp lammps-10Apr08/src/velocity.cpp
--- lammps-9Apr08/src/velocity.cpp	2008-01-18 10:18:46.000000000 -0700
+++ lammps-10Apr08/src/velocity.cpp	2008-04-08 16:55:40.000000000 -0600
@@ -197,35 +197,20 @@
       atom->map_set();
     }
 
-    random = new RanPark(lmp,seed);
+    // error check
 
     if (atom->tag_enable == 0)
       error->all("Cannot use velocity create loop all unless atoms have IDs");
-    int natoms = static_cast<int> (atom->natoms);
-
-    // check that atom IDs span range from 1 to natoms
+    if (atom->tag_consecutive() == 0)
+      error->all("Atom IDs must be consecutive for velocity create loop all");
 
-    int *tag = atom->tag;
-    int idmin = natoms;
-    int idmax = 0;
-
-    for (i = 0; i < nlocal; i++) {
-      idmin = MIN(idmin,tag[i]);
-      idmax = MAX(idmax,tag[i]);
-    }
-    int idminall,idmaxall;
-    MPI_Allreduce(&idmin,&idminall,1,MPI_INT,MPI_MIN,world);
-    MPI_Allreduce(&idmax,&idmaxall,1,MPI_INT,MPI_MAX,world);
-
-    if (idminall != 1 || idmaxall != natoms) {
-      char *str = (char *) "Cannot use velocity create loop all with non-contiguous atom IDs";
-      error->all(str);
-    }
-    
     // loop over all atoms in system
     // generate RNGs for all atoms, only assign to ones I own
     // use either per-type mass or per-atom rmass
 
+    random = new RanPark(lmp,seed);
+    int natoms = static_cast<int> (atom->natoms);
+
     for (i = 1; i <= natoms; i++) {
       if (dist_flag == 0) {
 	vx = random->uniform();