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

src/ASPHERE/fix_npt_asphere.cpp

@@ -44,14 +44,6 @@ FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) :
 	       "quat, angmom, torque, shape");
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixNPTAsphere::init()
-{
-  FixNPT::init();
-  dtq = 0.5 * update->dt;
-}
-
 /* ----------------------------------------------------------------------
    1st half of Verlet update 
 ------------------------------------------------------------------------- */
@@ -61,6 +53,8 @@ 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;
 
@@ -249,14 +243,6 @@ void FixNPTAsphere::final_integrate()
   }
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixNPTAsphere::reset_dt()
-{
-  FixNPT::reset_dt();
-  dtq = 0.5 * update->dt;
-}
-
 /* ----------------------------------------------------------------------
    Richardson iteration to update quaternion accurately
 ------------------------------------------------------------------------- */

src/ASPHERE/fix_npt_asphere.h

@@ -21,10 +21,8 @@ namespace LAMMPS_NS {
 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;

src/ASPHERE/fix_nve_asphere.cpp

@@ -66,6 +66,8 @@ void FixNVEAsphere::initial_integrate(int vflag)
 {
   double dtfm;
 
+  dtq = 0.5 * dtv;
+
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;

src/ASPHERE/fix_nve_asphere.h

@@ -27,6 +27,7 @@ class FixNVEAsphere : public FixNVE {
   void final_integrate();
 
  private:
+  double dtq;
   double **inertia;
 
   void richardson(double *, double *, double *);

src/ASPHERE/fix_nvt_asphere.cpp

@@ -47,18 +47,12 @@ FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-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 @@ void FixNVTAsphere::final_integrate()
   eta_dot *= drag_factor;
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixNVTAsphere::reset_dt()
-{
-  FixNVT::reset_dt();
-  dtq = 0.5 * update->dt;
-}
-
 /* ----------------------------------------------------------------------
    Richardson iteration to update quaternion accurately
 ------------------------------------------------------------------------- */

src/ASPHERE/fix_nvt_asphere.h

@@ -21,10 +21,8 @@ namespace LAMMPS_NS {
 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;

src/atom.cpp

@@ -601,26 +601,29 @@ void Atom::tag_extend()
 }
 
 /* ----------------------------------------------------------------------
-   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;
-  for (int i = 0; i < nlocal; i++) {
-    if (tag[i] <= 0) check[0] = 1;
-    if (tag[i] > check[1]) check[1] = tag[i];
-  }
-  MPI_Allreduce(check,check_all,2,MPI_INT,MPI_MAX,world);
 
-  if (check_all[0] || check_all[1] != static_cast<int> (natoms)) return 0;
+  int idmin = static_cast<int> (natoms);
+  int idmax = 0;
+  
+  for (int 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 != static_cast<int> (natoms)) return 0;
   return 1;
 }
 

src/fix_npt.h

@@ -25,7 +25,7 @@ class FixNPT : public Fix {
   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 @@ class FixNPT : public Fix {
   void write_restart(FILE *);
   void restart(char *);
   int modify_param(int, char **);
-  virtual void reset_dt();
+  void reset_dt();
 
  protected:
   int dimension,which;

src/fix_npt_sphere.cpp

@@ -53,18 +53,14 @@ FixNPTSphere::~FixNPTSphere()
 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 @@ void FixNPTSphere::initial_integrate(int vflag)
     }
   }
 
+  // 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 @@ void FixNPTSphere::final_integrate()
   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 @@ void FixNPTSphere::final_integrate()
     omega_dot[i] *= drag_factor;
   }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixNPTSphere::reset_dt()
-{
-  FixNPT::reset_dt();
-  dtfrotate = dtf / INERTIA;
-}

src/fix_npt_sphere.h

@@ -25,10 +25,8 @@ class FixNPTSphere : public FixNPT {
   void init();
   void initial_integrate(int);
   void final_integrate();
-  void reset_dt();
 
  private:
-  double dtfrotate;
   double *dttype;
   double factor_rotate;
 };

src/fix_nve.cpp

@@ -50,7 +50,6 @@ void FixNVE::init()
 {
   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 @@ void FixNVE::reset_dt()
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
 }

src/fix_nve.h

@@ -30,7 +30,7 @@ class FixNVE : public Fix {
   void reset_dt();
 
  protected:
-  double dtv,dtf,dtq;
+  double dtv,dtf;
   double *step_respa;
   int mass_require;
 };

src/fix_nve_sphere.cpp

@@ -31,7 +31,7 @@ enum{NONE,DIPOLE};
 /* ---------------------------------------------------------------------- */
 
 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 @@ void FixNVESphere::init()
 {
   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 @@ void FixNVESphere::init()
     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 @@ void FixNVESphere::initial_integrate(int vflag)
   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 @@ void FixNVESphere::final_integrate()
   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 @@ void FixNVESphere::final_integrate()
     }
     
   } 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::final_integrate()
     }
   }
 }
-
-/* ---------------------------------------------------------------------- */
-
-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;
-}

src/fix_nve_sphere.h

@@ -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 @@ class FixNVESphere : public Fix {
   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;
 };
 

src/fix_nvt.h

@@ -34,7 +34,7 @@ class FixNVT : public Fix {
   void restart(char *);
   int modify_param(int, char **);
   void reset_target(double);
-  virtual void reset_dt();
+  void reset_dt();
 
  protected:
   int which;

src/fix_nvt_sphere.cpp

@@ -54,18 +54,14 @@ FixNVTSphere::~FixNVTSphere()
 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 @@ void FixNVTSphere::initial_integrate(int vflag)
   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 @@ void FixNVTSphere::final_integrate()
   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 @@ void FixNVTSphere::final_integrate()
   eta_dot += f_eta*dthalf;
   eta_dot *= drag_factor;
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixNVTSphere::reset_dt()
-{
-  FixNVT::reset_dt();
-  dtfrotate = dtf / INERTIA;
-}

src/fix_nvt_sphere.h

@@ -25,10 +25,8 @@ class FixNVTSphere : public FixNVT {
   void init();
   void initial_integrate(int);
   void final_integrate();
-  void reset_dt();
 
  private:
-  double dtfrotate;
   double *dttype;
 };
 

src/velocity.cpp

@@ -197,35 +197,20 @@ void Velocity::create(int narg, char **arg)
       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);
+    if (atom->tag_consecutive() == 0)
+      error->all("Atom IDs must be consecutive for velocity create loop all");
 
-    // check that atom IDs span range from 1 to natoms
-
-    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();