Merge pull request #1002 from rbberger/set_velocity_command

Allow 'set' command to change atom velocities

doc/src/set.txt

@@ -36,6 +36,8 @@ keyword = {type} or {type/fraction} or {mol} or {x} or {y} or {z} or \
     value can be an atom-style variable (see below)
   {x},{y},{z} value = atom coordinate (distance units)
     value can be an atom-style variable (see below)
+  {vx},{vy},{vz} value = atom velocity (velocity units)
+    value can be an atom-style variable (see below)
   {charge} value = atomic charge (charge units)
     value can be an atom-style variable (see below)
   {dipole} values = x y z
@@ -127,6 +129,7 @@ set type 3 charge 0.5
 set type 1*3 charge 0.5
 set atom * charge v_atomfile
 set atom 100*200 x 0.5 y 1.0
+set atom 100 vx 0.0 vy 0.0 vz -1.0
 set atom 1492 type 3 :pre
 
 [Description:]
@@ -225,7 +228,8 @@ IDs.
 
 Keywords {x}, {y}, {z}, and {charge} set the coordinates or charge of
 all selected atoms.  For {charge}, the "atom style"_atom_style.html
-being used must support the use of atomic charge.
+being used must support the use of atomic charge. Keywords {vx}, {vy},
+and {vz} set the velocities of all selected atoms.
 
 Keyword {dipole} uses the specified x,y,z values as components of a
 vector to set as the orientation of the dipole moment vectors of the

python/lammps.py

@@ -708,6 +708,12 @@ class Atom(object):
             self.lmp.eval("vy[%d]" % self.index),
             self.lmp.eval("vz[%d]" % self.index))
 
+  @velocity.setter
+  def velocity(self, value):
+     self.lmp.set("atom", self.index, "vx", value[0])
+     self.lmp.set("atom", self.index, "vy", value[1])
+     self.lmp.set("atom", self.index, "vz", value[2])
+
   @property
   def force(self):
     return (self.lmp.eval("fx[%d]" % self.index),
@@ -738,6 +744,11 @@ class Atom2D(Atom):
     return (self.lmp.eval("vx[%d]" % self.index),
             self.lmp.eval("vy[%d]" % self.index))
 
+  @velocity.setter
+  def velocity(self, value):
+     self.lmp.set("atom", self.index, "vx", value[0])
+     self.lmp.set("atom", self.index, "vy", value[1])
+
   @property
   def force(self):
     return (self.lmp.eval("fx[%d]" % self.index),

src/set.cpp

@@ -48,7 +48,7 @@ enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
      THETA,THETA_RANDOM,ANGMOM,OMEGA,
      DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
      MESO_E,MESO_CV,MESO_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
-     SMD_CONTACT_RADIUS,DPDTHETA,INAME,DNAME};
+     SMD_CONTACT_RADIUS,DPDTHETA,INAME,DNAME,VX,VY,VZ};
 
 #define BIG INT_MAX
 
@@ -141,6 +141,27 @@ void Set::command(int narg, char **arg)
       set(Z);
       iarg += 2;
 
+    } else if (strcmp(arg[iarg],"vx") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
+      else dvalue = force->numeric(FLERR,arg[iarg+1]);
+      set(VX);
+      iarg += 2;
+
+    } else if (strcmp(arg[iarg],"vy") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
+      else dvalue = force->numeric(FLERR,arg[iarg+1]);
+      set(VY);
+      iarg += 2;
+
+    } else if (strcmp(arg[iarg],"vz") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
+      else dvalue = force->numeric(FLERR,arg[iarg+1]);
+      set(VZ);
+      iarg += 2;
+
     } else if (strcmp(arg[iarg],"charge") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
@@ -732,6 +753,9 @@ void Set::set(int keyword)
     else if (keyword == X) atom->x[i][0] = dvalue;
     else if (keyword == Y) atom->x[i][1] = dvalue;
     else if (keyword == Z) atom->x[i][2] = dvalue;
+    else if (keyword == VX) atom->v[i][0] = dvalue;
+    else if (keyword == VY) atom->v[i][1] = dvalue;
+    else if (keyword == VZ) atom->v[i][2] = dvalue;
     else if (keyword == CHARGE) atom->q[i] = dvalue;
     else if (keyword == MASS) {
       if (dvalue <= 0.0) error->one(FLERR,"Invalid mass in set command");