From 3b2784bc163fb9e467a5cd558aa1b2e71fa2ad8f Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 10 Jan 2011 17:07:41 -0500
Subject: [PATCH] implement a "reverse" sub-command for velocity.

---
 doc/velocity.txt |  8 +++++++-
 src/velocity.cpp | 22 +++++++++++++++++++++-
 src/velocity.h   |  1 +
 3 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/doc/velocity.txt b/doc/velocity.txt
index 0fbd0da8d1..1dee3954b5 100644
--- a/doc/velocity.txt
+++ b/doc/velocity.txt
@@ -13,7 +13,7 @@ velocity command :h3
 velocity group-ID style args keyword value ... :pre
 
 group-ID = ID of group of atoms whose velocity will be changed :ulb,l
-style = {create} or {set} or {scale} or {ramp} or {zero} :l
+style = {create} or {set} or {scale} or {ramp} or {reverse} or {zero} :l
   {create} args = temp seed
     temp = temperature value (temperature units)
     seed = random # seed (positive integer)
@@ -26,6 +26,7 @@ style = {create} or {set} or {scale} or {ramp} or {zero} :l
     vlo,vhi = lower and upper velocity value (velocity units)
     dim = {x} or {y} or {z}
     clo,chi = lower and upper coordinate bound (distance units)
+  {reverse}
   {zero} arg = {linear} or {angular}
     {linear} = zero the linear momentum
     {angular} = zero the angular momentum :pre
@@ -47,6 +48,7 @@ velocity all create 300.0 4928459 rot yes dist gaussian
 velocity border set NULL 4.0 3.0 sum yes units box
 velocity flow scale 300.0
 velocity flow ramp vx 0.0 5.0 y 5 25 temp mytemp
+velocity all reverse
 velocity all zero linear :pre
 
 [Description:]
@@ -77,6 +79,10 @@ outside the coordinate bounds (less than 5 or greater than 25 in this
 case), are assigned velocities equal to vlo or vhi (0.0 or 5.0 in this
 case).
 
+The {reverse} style replaces the velocities of the group of atoms with
+the negative of the original values. The total kinetic energy is unchanged.
+No other changes are made to the velocities of the atoms.
+
 The {zero} style adjusts the velocities of the group of atoms so that
 the aggregate linear or angular momentum is zero.  No other changes
 are made to the velocities of the atoms.
diff --git a/src/velocity.cpp b/src/velocity.cpp
index 4cadf22de6..12f44ce3c7 100644
--- a/src/velocity.cpp
+++ b/src/velocity.cpp
@@ -34,7 +34,7 @@
 
 using namespace LAMMPS_NS;
 
-enum{CREATE,SET,SCALE,RAMP,ZERO};
+enum{CREATE,SET,SCALE,RAMP,REVERSE,ZERO};
 enum{ALL,LOCAL,GEOM};
 
 #define WARMUP 100
@@ -74,6 +74,7 @@ void Velocity::command(int narg, char **arg)
   else if (strcmp(arg[1],"set") == 0) style = SET;
   else if (strcmp(arg[1],"scale") == 0) style = SCALE;
   else if (strcmp(arg[1],"ramp") == 0) style = RAMP;
+  else if (strcmp(arg[1],"reverse") == 0) style = REVERSE;
   else if (strcmp(arg[1],"zero") == 0) style = ZERO;
   else error->all("Illegal velocity command");
 
@@ -121,6 +122,7 @@ void Velocity::command(int narg, char **arg)
   else if (style == SET) set(narg-2,&arg[2]);
   else if (style == SCALE) scale(narg-2,&arg[2]);
   else if (style == RAMP) ramp(narg-2,&arg[2]);
+  else if (style == REVERSE) reverse(narg-2,&arg[2]);
   else if (style == ZERO) zero(narg-2,&arg[2]);
 }
 
@@ -528,6 +530,24 @@ void Velocity::rescale(double t_old, double t_new)
     }
 }
 
+/* ----------------------------------------------------------------------
+   reverse velocities, so an MD can run backwards.
+------------------------------------------------------------------------- */
+
+void Velocity::reverse(int narg, char **arg)
+{
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      v[i][0] = -v[i][0];
+      v[i][1] = -v[i][1];
+      v[i][2] = -v[i][2];
+    }
+}
+
 /* ----------------------------------------------------------------------
    zero the linear momentum of a group of atoms by adjusting v by -Vcm
 ------------------------------------------------------------------------- */
diff --git a/src/velocity.h b/src/velocity.h
index 3401663c2d..1850121293 100644
--- a/src/velocity.h
+++ b/src/velocity.h
@@ -39,6 +39,7 @@ class Velocity : protected Pointers {
   double xscale,yscale,zscale;
   class Compute *temperature;
 
+  void reverse(int, char **);
   void set(int, char **);
   void scale(int, char **);
   void ramp(int, char **);