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

src/compute.h

@@ -114,6 +114,7 @@ class Compute : protected Pointers {
   virtual int dof_remove(int) {return 0;}
   virtual void remove_bias(int, double *) {}
   virtual void remove_bias_all() {}
+  virtual void reapply_bias_all() {}
   virtual void restore_bias(int, double *) {}
   virtual void restore_bias_all() {}
 

src/compute_temp_partial.cpp

@@ -217,6 +217,38 @@ void ComputeTempPartial::remove_bias_all()
   }
 }
 
+/* ----------------------------------------------------------------------
+   reset thermal velocity of atoms to be consistent with bias
+   called from velocity command after creating thermal velocities
+   needed to re-zero components that should stay zero
+------------------------------------------------------------------------- */
+
+void ComputeTempPartial::reapply_bias_all()
+{
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (nlocal > maxbias) {
+    memory->destroy(vbiasall);
+    maxbias = atom->nmax;
+    memory->create(vbiasall,maxbias,3,"temp/partial:vbiasall");
+  }
+
+  if (!xflag) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) v[i][0] = 0.0;
+  }
+  if (!yflag) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) v[i][1] = 0.0;
+  }
+  if (!zflag) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) v[i][2] = 0.0;
+  }
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to atom I removed by remove_bias()
    assume remove_bias() was previously called

src/compute_temp_partial.h

@@ -36,6 +36,7 @@ class ComputeTempPartial : public Compute {
   int dof_remove(int);
   void remove_bias(int, double *);
   void remove_bias_all();
+  void reapply_bias_all();
   void restore_bias(int, double *);
   void restore_bias_all();
   double memory_usage();

src/velocity.cpp

@@ -84,6 +84,7 @@ void Velocity::command(int narg, char **arg)
   sum_flag = 0;
   momentum_flag = 1;
   rotation_flag = 0;
+  bias_flag = 0;
   loop_flag = ALL;
   scale_flag = 1;
   rfix = -1;
@@ -135,41 +136,59 @@ void Velocity::init_external(const char *extgroup)
 void Velocity::create(double t_desired, int seed)
 {
   int i;
+  double **vhold;
 
   if (seed <= 0) error->all(FLERR,"Illegal velocity create command");
 
-  // if temperature = NULL, create a new ComputeTemp with the velocity group
+  // if sum_flag set, store a copy of current velocities
 
-  int tflag = 0;
+  if (sum_flag) {
-  if (temperature == NULL) {
+    double **v = atom->v;
+    int nlocal = atom->nlocal;
+    memory->create(vhold,nlocal,3,"velocity:vhold");
+    for (i = 0; i < nlocal; i++) {
+      vhold[i][0] = v[i][0];
+      vhold[i][1] = v[i][1];
+      vhold[i][2] = v[i][2];
+    }
+  }
+
+  // if temperature = NULL or bias_flag set,
+  // create a new ComputeTemp with the velocity group
+
+  int tcreate_flag = 0;
+  Compute *temperature_nobias = NULL;
+
+  if (temperature == NULL || bias_flag) {
     char **arg = new char*[3];
     arg[0] = (char *) "velocity_temp";
     arg[1] = group->names[igroup];
     arg[2] = (char *) "temp";
-    temperature = new ComputeTemp(lmp,3,arg);
+    if (temperature == NULL) {
-    tflag = 1;
+      temperature = new ComputeTemp(lmp,3,arg);
+      tcreate_flag = 1;
+    } else temperature_nobias = new ComputeTemp(lmp,3,arg);
     delete [] arg;
   }
 
-  // initialize temperature computation
+  // initialize temperature computation(s)
   // warn if groups don't match
 
   if (igroup != temperature->igroup && comm->me == 0)
     error->warning(FLERR,"Mismatch between velocity and compute groups");
   temperature->init();
   temperature->setup();
+  if (temperature_nobias) {
+    temperature_nobias->init();
+    temperature_nobias->setup();
+  }
 
-  // store a copy of current velocities
+  // if bias_flag set, remove bias velocity from all atoms
+  // for some temperature computes, must first calculate temp to do that
 
-  double **v = atom->v;
+  if (bias_flag) {
-  int nlocal = atom->nlocal;
+    temperature->compute_scalar();
-  double **vhold;
+    temperature->remove_bias_all();
-  memory->create(vhold,nlocal,3,"velocity:vnew");
-
-  for (i = 0; i < nlocal; i++) {
-    vhold[i][0] = v[i][0];
-    vhold[i][1] = v[i][1];
-    vhold[i][2] = v[i][2];
   }
 
   // create new velocities, in uniform or gaussian distribution
@@ -185,13 +204,18 @@ void Velocity::create(double t_desired, int seed)
   //      via random->reset()
   //    will always produce same V, independent of P
   // adjust by factor for atom mass
-  // for 2d, set Vz to 0.0
+  // set xdim,ydim,zdim = 1/0 for whether to create velocity in those dims
+  //   zdim = 0 for 2d
+  //   any dims can be 0 if bias temperature compute turns them off
+  //     currently only temp/partial does
 
+  double **v = atom->v;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
-  int dimension = domain->dimension;
+  int nlocal = atom->nlocal;
+  int dim = domain->dimension;
 
   int m;
   double vx,vy,vz,factor;
@@ -244,7 +268,7 @@ void Velocity::create(double t_desired, int seed)
           else factor = 1.0/sqrt(mass[type[m]]);
           v[m][0] = vx * factor;
           v[m][1] = vy * factor;
-          if (dimension == 3) v[m][2] = vz * factor;
+          if (dim == 3) v[m][2] = vy * factor;
           else v[m][2] = 0.0;
         }
       }
@@ -276,7 +300,7 @@ void Velocity::create(double t_desired, int seed)
         else factor = 1.0/sqrt(mass[type[i]]);
         v[i][0] = vx * factor;
         v[i][1] = vy * factor;
-        if (dimension == 3) v[i][2] = vz * factor;
+        if (dim == 3) v[i][2] = vz * factor;
         else v[i][2] = 0.0;
       }
     }
@@ -302,7 +326,7 @@ void Velocity::create(double t_desired, int seed)
         else factor = 1.0/sqrt(mass[type[i]]);
         v[i][0] = vx * factor;
         v[i][1] = vy * factor;
-        if (dimension == 3) v[i][2] = vz * factor;
+        if (dim == 3) v[i][2] = vz * factor;
         else v[i][2] = 0.0;
       }
     }
@@ -314,10 +338,22 @@ void Velocity::create(double t_desired, int seed)
   if (rotation_flag) zero_rotation();
 
   // scale temp to desired value
+  // if bias flag is set, bias velocities have already been removed:
+  //   no-bias compute calculates temp only for new thermal velocities
 
-  double t = temperature->compute_scalar();
+  double t;
+  if (bias_flag == 0) t = temperature->compute_scalar();
+  else t = temperature_nobias->compute_scalar();
   rescale(t,t_desired);
 
+  // if bias_flag set, restore bias velocity to all atoms
+  // reapply for bias = compute temp/partial to reset v dims to 0.0
+
+  if (bias_flag) {
+    temperature->reapply_bias_all();
+    temperature->restore_bias_all();
+  }
+
   // if sum_flag set, add back in previous velocities
 
   if (sum_flag) {
@@ -328,14 +364,15 @@ void Velocity::create(double t_desired, int seed)
         v[i][2] += vhold[i][2];
       }
     }
+    memory->destroy(vhold);
   }
 
   // free local memory
-  // if temperature was created, delete it
+  // if temperature compute was created, delete it
 
-  memory->destroy(vhold);
   delete random;
-  if (tflag) delete temperature;
+  if (tcreate_flag) delete temperature;
+  if (temperature_nobias) delete temperature_nobias;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -535,9 +572,19 @@ void Velocity::scale(int narg, char **arg)
   temperature->setup();
 
   // scale temp to desired value
+  // if bias flag is set:
+  //   temperature calculation will be done accounting for bias
+  //   remove/restore bias velocities before/after rescale
 
-  double t = temperature->compute_scalar();
+  if (bias_flag == 0) {
-  rescale(t,t_desired);
+    double t = temperature->compute_scalar();
+    rescale(t,t_desired);
+  } else {
+    double t = temperature->compute_scalar();
+    temperature->remove_bias_all();
+    rescale(t,t_desired);
+    temperature->restore_bias_all();
+  }
 
   // if temperature was created, delete it
 
@@ -660,6 +707,7 @@ void Velocity::zero(int narg, char **arg)
 
 /* ----------------------------------------------------------------------
    rescale velocities of group atoms to t_new from t_old
+   no bias applied here, since done in create() and scale()
 ------------------------------------------------------------------------- */
 
 void Velocity::rescale(double t_old, double t_new)
@@ -804,6 +852,12 @@ void Velocity::options(int narg, char **arg)
         error->all(FLERR,
                    "Velocity temperature ID does not compute temperature");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"bias") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
+      if (strcmp(arg[iarg+1],"no") == 0) bias_flag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) bias_flag = 1;
+      else error->all(FLERR,"Illegal velocity command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"loop") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
       if (strcmp(arg[iarg+1],"all") == 0) loop_flag = ALL;
@@ -824,4 +878,11 @@ void Velocity::options(int narg, char **arg)
       iarg += 2;
     } else error->all(FLERR,"Illegal velocity command");
   }
+
+  // error check
+
+  if (bias_flag && temperature == NULL)
+    error->all(FLERR,"Cannot use velocity bias command without temp keyword");
+  if (bias_flag && temperature->tempbias == 0)
+    error->all(FLERR,"Velocity temperature ID does calculate a velocity bias");
 }

src/velocity.h

@@ -35,7 +35,8 @@ class Velocity : protected Pointers {
  private:
   int igroup,groupbit;
   int style;
-  int dist_flag,sum_flag,momentum_flag,rotation_flag,loop_flag,scale_flag,rfix;
+  int dist_flag,sum_flag,momentum_flag,rotation_flag;
+  int bias_flag,loop_flag,scale_flag,rfix;
   double xscale,yscale,zscale;
   class Compute *temperature;