diff -Naur lammps-11Apr08/doc/fix_nph.html lammps-12Apr08/doc/fix_nph.html
--- lammps-11Apr08/doc/fix_nph.html	2008-04-04 15:15:50.000000000 -0600
+++ lammps-12Apr08/doc/fix_nph.html	2008-04-10 14:53:47.000000000 -0600
@@ -125,8 +125,10 @@
 oscillations after a few periods.
 </P>
 <P>For all pressure styles, the simulation box stays rectangular in
-shape.  Parinello-Rahman boundary conditions (tilted box) for this fix
-are not yet implemented in LAMMPS.
+shape.  Parinello-Rahman boundary condition for tilted boxes
+(triclinic symmetry) are supported by other LAMMPS commands (see <A HREF = "Section_howto.html#4_12">this
+section</A> of the manual), but not yet by this
+command.
 </P>
 <P>For all styles, the <I>Pdamp</I> parameter determines the time scale on
 which pressure is relaxed.  For example, a value of 1000.0 means to
diff -Naur lammps-11Apr08/doc/fix_nph.txt lammps-12Apr08/doc/fix_nph.txt
--- lammps-11Apr08/doc/fix_nph.txt	2008-04-04 15:15:50.000000000 -0600
+++ lammps-12Apr08/doc/fix_nph.txt	2008-04-10 14:53:47.000000000 -0600
@@ -116,8 +116,10 @@
 oscillations after a few periods.
 
 For all pressure styles, the simulation box stays rectangular in
-shape.  Parinello-Rahman boundary conditions (tilted box) for this fix
-are not yet implemented in LAMMPS.
+shape.  Parinello-Rahman boundary condition for tilted boxes
+(triclinic symmetry) are supported by other LAMMPS commands (see "this
+section"_Section_howto.html#4_12 of the manual), but not yet by this
+command.
 
 For all styles, the {Pdamp} parameter determines the time scale on
 which pressure is relaxed.  For example, a value of 1000.0 means to
diff -Naur lammps-11Apr08/doc/fix_npt.html lammps-12Apr08/doc/fix_npt.html
--- lammps-11Apr08/doc/fix_npt.html	2008-04-04 15:15:50.000000000 -0600
+++ lammps-12Apr08/doc/fix_npt.html	2008-04-10 14:53:47.000000000 -0600
@@ -136,8 +136,10 @@
 damp oscillations after a few periods.
 </P>
 <P>For all pressure styles, the simulation box stays rectangular in
-shape.  Parinello-Rahman boundary conditions (tilted box) for this fix
-are not yet implemented in LAMMPS.
+shape.  Parinello-Rahman boundary condition for tilted boxes
+(triclinic symmetry) are supported by other LAMMPS commands (see <A HREF = "Section_howto.html#4_12">this
+section</A> of the manual), but not yet by this
+command.
 </P>
 <P>For all styles, the <I>Pdamp</I> parameter operates like the <I>Tdamp</I>
 parameter, determining the time scale on which pressure is relaxed.
diff -Naur lammps-11Apr08/doc/fix_npt.txt lammps-12Apr08/doc/fix_npt.txt
--- lammps-11Apr08/doc/fix_npt.txt	2008-04-04 15:15:50.000000000 -0600
+++ lammps-12Apr08/doc/fix_npt.txt	2008-04-10 14:53:47.000000000 -0600
@@ -125,8 +125,10 @@
 damp oscillations after a few periods.
 
 For all pressure styles, the simulation box stays rectangular in
-shape.  Parinello-Rahman boundary conditions (tilted box) for this fix
-are not yet implemented in LAMMPS.
+shape.  Parinello-Rahman boundary condition for tilted boxes
+(triclinic symmetry) are supported by other LAMMPS commands (see "this
+section"_Section_howto.html#4_12 of the manual), but not yet by this
+command.
 
 For all styles, the {Pdamp} parameter operates like the {Tdamp}
 parameter, determining the time scale on which pressure is relaxed.
diff -Naur lammps-11Apr08/doc/fix_press_berendsen.html lammps-12Apr08/doc/fix_press_berendsen.html
--- lammps-11Apr08/doc/fix_press_berendsen.html	2008-04-04 15:15:50.000000000 -0600
+++ lammps-12Apr08/doc/fix_press_berendsen.html	2008-04-10 14:53:47.000000000 -0600
@@ -117,8 +117,10 @@
 oscillations after a few periods.
 </P>
 <P>For all pressure styles, the simulation box stays rectangular in
-shape.  Parinello-Rahman boundary conditions (tilted box) for this fix
-are not yet implemented in LAMMPS.
+shape.  Parinello-Rahman boundary condition for tilted boxes
+(triclinic symmetry) are supported by other LAMMPS commands (see <A HREF = "Section_howto.html#4_12">this
+section</A> of the manual), but not yet by this
+command.
 </P>
 <P>For all styles, the <I>Pdamp</I> parameter determines the time scale on
 which pressure is relaxed.  For example, a value of 1000.0 means to
diff -Naur lammps-11Apr08/doc/fix_press_berendsen.txt lammps-12Apr08/doc/fix_press_berendsen.txt
--- lammps-11Apr08/doc/fix_press_berendsen.txt	2008-04-04 15:15:50.000000000 -0600
+++ lammps-12Apr08/doc/fix_press_berendsen.txt	2008-04-10 14:53:47.000000000 -0600
@@ -108,8 +108,10 @@
 oscillations after a few periods.
 
 For all pressure styles, the simulation box stays rectangular in
-shape.  Parinello-Rahman boundary conditions (tilted box) for this fix
-are not yet implemented in LAMMPS.
+shape.  Parinello-Rahman boundary condition for tilted boxes
+(triclinic symmetry) are supported by other LAMMPS commands (see "this
+section"_Section_howto.html#4_12 of the manual), but not yet by this
+command.
 
 For all styles, the {Pdamp} parameter determines the time scale on
 which pressure is relaxed.  For example, a value of 1000.0 means to
diff -Naur lammps-11Apr08/doc/units.html lammps-12Apr08/doc/units.html
--- lammps-11Apr08/doc/units.html	2008-02-29 18:13:20.000000000 -0700
+++ lammps-12Apr08/doc/units.html	2008-04-10 14:53:47.000000000 -0600
@@ -34,44 +34,51 @@
 www.physics.nist.gov. For the definition of Kcal in real units, LAMMPS
 uses the thermochemical calorie = 4.184 J.
 </P>
-<P>For style <I>lj</I>, all quantities are unitless:
-</P>
-<UL><LI>distance = sigma
-<LI>time = tau
-<LI>mass = one
-<LI>energy = epsilon
-<LI>velocity = sigma/tau
-<LI>force = epsilon/sigma
-<LI>temperature = reduced LJ temperature
-<LI>pressure = reduced LJ pressure
-<LI>charge = reduced LJ charge
-<LI>dipole = reduced LJ dipole moment
-<LI>electric field = force/charge 
+<P>For style <I>lj</I>, all quantities are unitless.  The formula relating the
+reduced or unitless quantity (with an asterisk) to the same quantity
+with units is also given:
+</P>
+<UL><LI>m (mass) = epsilon = sigma = tau = Boltzmann constant = 1
+<LI>distance = sigma, where x* = x / sigma
+<LI>time = tau, where tau = t* = t (Kb T / m / sigma^2)^1/2
+<LI>energy = epsilon, where E* = E / epsilon
+<LI>velocity = sigma/tau, where v* = v tau / sigma
+<LI>force = epsilon/sigma, where f* = f sigma / epsilon
+<LI>temperature = reduced LJ temperature, where T* = T Kb / epsilon
+<LI>pressure = reduced LJ pressure, where P* = P sigma^3 / epsilon
+<LI>viscosity = reduced LJ viscosity, where eta* = eta sigma^3 / epsilon / tau
+<LI>charge = reduced LJ charge, where q* = q / (4 pi perm0 sigma epsilon)^1/2
+<LI>dipole = reduced LJ dipole, moment where *mu = mu / (4 pi perm0 sigma^3 epsilon)^1/2
+<LI>electric field = force/charge, where E* = E (4 pi perm0 sigma epsilon)^1/2 sigma / epsilon 
 </UL>
 <P>For style <I>real</I>, these are the units:
 </P>
-<UL><LI>distance = Angstroms
-<LI>time = femtoseconds
+<UL><LI>Boltzmann constant = 0.0019872067 Kcal/mole per degree K
 <LI>mass = grams/mole
+<LI>distance = Angstroms
+<LI>time = femtoseconds
 <LI>energy = Kcal/mole 
 <LI>velocity = Angstroms/femtosecond 
 <LI>force = Kcal/mole-Angstrom
 <LI>temperature = degrees K
 <LI>pressure = atmospheres
+<LI>viscosity = Poise
 <LI>charge = multiple of electron charge (+1.0 is a proton)
 <LI>dipole = charge*Angstroms
 <LI>electric field = volts/Angstrom 
 </UL>
 <P>For style <I>metal</I>, these are the units:
 </P>
-<UL><LI>distance = Angstroms
-<LI>time = picoseconds
+<UL><LI>Boltzmann constant = 8.617343e-5 eV per degree K
 <LI>mass = grams/mole
+<LI>distance = Angstroms
+<LI>time = picoseconds
 <LI>energy = eV
 <LI>velocity = Angstroms/picosecond 
 <LI>force = eV/Angstrom
 <LI>temperature = degrees K
 <LI>pressure = bars
+<LI>viscosity = Poise
 <LI>charge = multiple of electron charge (+1.0 is a proton)
 <LI>dipole = charge*Angstroms
 <LI>electric field = volts/Angstrom 
diff -Naur lammps-11Apr08/doc/units.txt lammps-12Apr08/doc/units.txt
--- lammps-11Apr08/doc/units.txt	2008-02-29 18:13:20.000000000 -0700
+++ lammps-12Apr08/doc/units.txt	2008-04-10 14:53:47.000000000 -0600
@@ -31,44 +31,51 @@
 www.physics.nist.gov. For the definition of Kcal in real units, LAMMPS
 uses the thermochemical calorie = 4.184 J.
 
-For style {lj}, all quantities are unitless:
-
-distance = sigma
-time = tau
-mass = one
-energy = epsilon
-velocity = sigma/tau
-force = epsilon/sigma
-temperature = reduced LJ temperature
-pressure = reduced LJ pressure
-charge = reduced LJ charge
-dipole = reduced LJ dipole moment
-electric field = force/charge :ul
+For style {lj}, all quantities are unitless.  The formula relating the
+reduced or unitless quantity (with an asterisk) to the same quantity
+with units is also given:
+
+m (mass) = epsilon = sigma = tau = Boltzmann constant = 1
+distance = sigma, where x* = x / sigma
+time = tau, where tau = t* = t (Kb T / m / sigma^2)^1/2
+energy = epsilon, where E* = E / epsilon
+velocity = sigma/tau, where v* = v tau / sigma
+force = epsilon/sigma, where f* = f sigma / epsilon
+temperature = reduced LJ temperature, where T* = T Kb / epsilon
+pressure = reduced LJ pressure, where P* = P sigma^3 / epsilon
+viscosity = reduced LJ viscosity, where eta* = eta sigma^3 / epsilon / tau
+charge = reduced LJ charge, where q* = q / (4 pi perm0 sigma epsilon)^1/2
+dipole = reduced LJ dipole, moment where *mu = mu / (4 pi perm0 sigma^3 epsilon)^1/2
+electric field = force/charge, where E* = E (4 pi perm0 sigma epsilon)^1/2 sigma / epsilon :ul
 
 For style {real}, these are the units:
 
+Boltzmann constant = 0.0019872067 Kcal/mole per degree K
+mass = grams/mole
 distance = Angstroms
 time = femtoseconds
-mass = grams/mole
 energy = Kcal/mole 
 velocity = Angstroms/femtosecond 
 force = Kcal/mole-Angstrom
 temperature = degrees K
 pressure = atmospheres
+viscosity = Poise
 charge = multiple of electron charge (+1.0 is a proton)
 dipole = charge*Angstroms
 electric field = volts/Angstrom :ul
 
 For style {metal}, these are the units:
 
+Boltzmann constant = 8.617343e-5 eV per degree K
+mass = grams/mole
 distance = Angstroms
 time = picoseconds
-mass = grams/mole
 energy = eV
 velocity = Angstroms/picosecond 
 force = eV/Angstrom
 temperature = degrees K
 pressure = bars
+viscosity = Poise
 charge = multiple of electron charge (+1.0 is a proton)
 dipole = charge*Angstroms
 electric field = volts/Angstrom :ul
diff -Naur lammps-11Apr08/src/COLLOID/pair_lubricate.cpp lammps-12Apr08/src/COLLOID/pair_lubricate.cpp
--- lammps-11Apr08/src/COLLOID/pair_lubricate.cpp	2008-02-14 16:05:04.000000000 -0700
+++ lammps-12Apr08/src/COLLOID/pair_lubricate.cpp	2008-04-10 14:53:37.000000000 -0600
@@ -86,6 +86,7 @@
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
   int omega_flag = atom->omega_flag;
+  double vxmu2f = force->vxmu2f;
 
   inum = list->inum;
   ilist = list->ilist;
@@ -198,16 +199,17 @@
         h_sep = r - 2.0*radi;
 
 	if (flag1)
-	  a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); 
+	  a_squeeze = vxmu2f * 
+	    (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); 
 	if (flag2) 
-	  a_shear = (PI*mu*2.*radi/2.0) *
+	  a_shear = vxmu2f * (PI*mu*2.*radi/2.0) *
 	    log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0;
 	if (flag3) 
-	  a_pump = (PI*mu*pow(2.0*radi,4)/8.0) *
+	  a_pump = vxmu2f * (PI*mu*pow(2.0*radi,4)/8.0) *
 	    ((3.0/20.0) * log(2.0*radi/2.0/h_sep) + 
 	     (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep));
 	if (flag4)
-	  a_twist = (PI*mu*pow(2.0*radi,4)/4.0) *
+	  a_twist = vxmu2f * (PI*mu*pow(2.0*radi,4)/4.0) *
 	    (h_sep/2.0/radi) * log(2.0/(2.0*h_sep));
 
         if (h_sep >= cut_inner[itype][jtype]) {
diff -Naur lammps-11Apr08/src/fix_nph.cpp lammps-12Apr08/src/fix_nph.cpp
--- lammps-11Apr08/src/fix_nph.cpp	2008-04-04 17:18:14.000000000 -0600
+++ lammps-12Apr08/src/fix_nph.cpp	2008-04-10 14:53:37.000000000 -0600
@@ -37,6 +37,8 @@
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
 
+enum{XYZ,XY,YZ,XZ,ANISO};
+
 /* ---------------------------------------------------------------------- */
 
 FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
@@ -55,7 +57,7 @@
   if (strcmp(arg[3],"xyz") == 0) {
     if (narg < 7) error->all("Illegal fix nph command");
 
-    press_couple = 0;
+    press_couple = XYZ;
     p_start[0] = p_start[1] = p_start[2] = atof(arg[4]);
     p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]);
     p_period[0] = p_period[1] = p_period[2] = atof(arg[6]);
@@ -66,16 +68,16 @@
     }
 
   } else {
-    if (strcmp(arg[3],"xy") == 0) press_couple = 1;
-    else if (strcmp(arg[3],"yz") == 0) press_couple = 2;
-    else if (strcmp(arg[3],"xz") == 0) press_couple = 3;
-    else if (strcmp(arg[3],"aniso") == 0) press_couple = 4;
+    if (strcmp(arg[3],"xy") == 0) press_couple = XY;
+    else if (strcmp(arg[3],"yz") == 0) press_couple = YZ;
+    else if (strcmp(arg[3],"xz") == 0) press_couple = XZ;
+    else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO;
     else error->all("Illegal fix nph command");
 
     if (narg < 11) error->all("Illegal fix nph command");
 
     if (domain->dimension == 2 && 
-	(press_couple == 1 || press_couple == 2 || press_couple == 3))
+	(press_couple == XY || press_couple == YZ || press_couple == XZ))
       error->all("Invalid fix nph command for a 2d simulation");
 
     if (strcmp(arg[4],"NULL") == 0) {
@@ -117,7 +119,7 @@
   allremap = 1;
 
   int iarg;
-  if (press_couple == 0) iarg = 7;
+  if (press_couple == XYZ) iarg = 7;
   else iarg = 11;
 
   while (iarg < narg) {
@@ -134,7 +136,7 @@
     } else error->all("Illegal fix nph command");
   }
 
-  // error checks
+  // check for periodicity in controlled dimensions
 
   if (p_flag[0] && domain->xperiodic == 0)
     error->all("Cannot use fix nph on a non-periodic dimension");
@@ -313,7 +315,7 @@
   p_target[1] = p_start[1];
   p_target[2] = p_start[2];
 
-  if (press_couple == 0) {
+  if (press_couple == XYZ) {
     double tmp = temperature->compute_scalar();
     tmp = pressure->compute_scalar();
   } else {
@@ -426,7 +428,7 @@
 
   // compute new pressure
 
-  if (press_couple == 0) {
+  if (press_couple == XYZ) {
     double tmp = temperature->compute_scalar();
     tmp = pressure->compute_scalar();
   } else {
@@ -598,21 +600,21 @@
 {
   double *tensor = pressure->vector;
 
-  if (press_couple == 0)
+  if (press_couple == XYZ)
     p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
-  else if (press_couple == 1) {
+  else if (press_couple == XY) {
     double ave = 0.5 * (tensor[0] + tensor[1]);
     p_current[0] = p_current[1] = ave;
     p_current[2] = tensor[2];
-  } else if (press_couple == 2) {
+  } else if (press_couple == YZ) {
     double ave = 0.5 * (tensor[1] + tensor[2]);
     p_current[1] = p_current[2] = ave;
     p_current[0] = tensor[0];
-  } else if (press_couple == 3) {
+  } else if (press_couple == XZ) {
     double ave = 0.5 * (tensor[0] + tensor[2]);
     p_current[0] = p_current[2] = ave;
     p_current[1] = tensor[1];
-  } if (press_couple == 4) {
+  } else if (press_couple == ANISO) {
     p_current[0] = tensor[0];
     p_current[1] = tensor[1];
     p_current[2] = tensor[2];
diff -Naur lammps-11Apr08/src/fix_npt.cpp lammps-12Apr08/src/fix_npt.cpp
--- lammps-11Apr08/src/fix_npt.cpp	2008-04-04 17:18:14.000000000 -0600
+++ lammps-12Apr08/src/fix_npt.cpp	2008-04-10 14:53:37.000000000 -0600
@@ -37,6 +37,7 @@
 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
 
 enum{NOBIAS,BIAS};
+enum{XYZ,XY,YZ,XZ,ANISO};
 
 /* ---------------------------------------------------------------------- */
 
@@ -63,7 +64,7 @@
   if (strcmp(arg[6],"xyz") == 0) {
     if (narg < 10) error->all("Illegal fix npt command");
 
-    press_couple = 0;
+    press_couple = XYZ;
     p_start[0] = p_start[1] = p_start[2] = atof(arg[7]);
     p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]);
     p_period[0] = p_period[1] = p_period[2] = atof(arg[9]);
@@ -74,16 +75,16 @@
     }
 
   } else {
-    if (strcmp(arg[6],"xy") == 0) press_couple = 1;
-    else if (strcmp(arg[6],"yz") == 0) press_couple = 2;
-    else if (strcmp(arg[6],"xz") == 0) press_couple = 3;
-    else if (strcmp(arg[6],"aniso") == 0) press_couple = 4;
+    if (strcmp(arg[6],"xy") == 0) press_couple = XY;
+    else if (strcmp(arg[6],"yz") == 0) press_couple = YZ;
+    else if (strcmp(arg[6],"xz") == 0) press_couple = XZ;
+    else if (strcmp(arg[6],"aniso") == 0) press_couple = ANISO;
     else error->all("Illegal fix npt command");
 
     if (narg < 14) error->all("Illegal fix npt command");
 
     if (domain->dimension == 2 && 
-	(press_couple == 1 || press_couple == 2 || press_couple == 3))
+	(press_couple == XY || press_couple == YZ || press_couple == XZ))
       error->all("Invalid fix npt command for a 2d simulation");
 
     if (strcmp(arg[7],"NULL") == 0) {
@@ -125,7 +126,7 @@
   allremap = 1;
 
   int iarg;
-  if (press_couple == 0) iarg = 10;
+  if (press_couple == XYZ) iarg = 10;
   else iarg = 14;
 
   while (iarg < narg) {
@@ -142,7 +143,7 @@
     } else error->all("Illegal fix npt command");
   }
 
-  // error checks
+  // check for periodicity in controlled dimensions
 
   if (p_flag[0] && domain->xperiodic == 0)
     error->all("Cannot use fix npt on a non-periodic dimension");
@@ -322,7 +323,7 @@
   p_target[2] = p_start[2];
 
   t_current = temperature->compute_scalar();
-  if (press_couple == 0) {
+  if (press_couple == XYZ) {
     double tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
@@ -473,7 +474,7 @@
   // compute new T,P
 
   t_current = temperature->compute_scalar();
-  if (press_couple == 0) {
+  if (press_couple == XYZ) {
     double tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
@@ -687,21 +688,21 @@
 {
   double *tensor = pressure->vector;
 
-  if (press_couple == 0)
+  if (press_couple == XYZ)
     p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
-  else if (press_couple == 1) {
+  else if (press_couple == XY) {
     double ave = 0.5 * (tensor[0] + tensor[1]);
     p_current[0] = p_current[1] = ave;
     p_current[2] = tensor[2];
-  } else if (press_couple == 2) {
+  } else if (press_couple == YZ) {
     double ave = 0.5 * (tensor[1] + tensor[2]);
     p_current[1] = p_current[2] = ave;
     p_current[0] = tensor[0];
-  } else if (press_couple == 3) {
+  } else if (press_couple == XZ) {
     double ave = 0.5 * (tensor[0] + tensor[2]);
     p_current[0] = p_current[2] = ave;
     p_current[1] = tensor[1];
-  } if (press_couple == 4) {
+  } else if (press_couple == ANISO) {
     p_current[0] = tensor[0];
     p_current[1] = tensor[1];
     p_current[2] = tensor[2];
diff -Naur lammps-11Apr08/src/force.h lammps-12Apr08/src/force.h
--- lammps-11Apr08/src/force.h	2007-10-04 11:57:04.000000000 -0600
+++ lammps-12Apr08/src/force.h	2008-04-10 14:53:37.000000000 -0600
@@ -26,6 +26,7 @@
   double nktv2p;                     // conversion of NkT/V to pressure
   double qqr2e;                      // conversion of q^2/r to energy
   double qe2f;                       // conversion of qE to force
+  double vxmu2f;                     // conversion of vx mu to force
   double dielectric;                 // dielectric constant
   double qqrd2e;                     // q^2/r to energy w/ dielectric constant
 
diff -Naur lammps-11Apr08/src/min_cg.cpp lammps-12Apr08/src/min_cg.cpp
--- lammps-11Apr08/src/min_cg.cpp	2008-01-09 14:56:57.000000000 -0700
+++ lammps-12Apr08/src/min_cg.cpp	2008-04-10 14:53:37.000000000 -0600
@@ -186,6 +186,7 @@
     eng_force(&ntmp,&xtmp,&htmp,&etmp);
     output->write(update->ntimestep);
   }
+
   timer->barrier_stop(TIME_LOOP);
 
   // delete fix at end of run, so its atom arrays won't persist
diff -Naur lammps-11Apr08/src/update.cpp lammps-12Apr08/src/update.cpp
--- lammps-11Apr08/src/update.cpp	2007-10-22 17:26:34.000000000 -0600
+++ lammps-12Apr08/src/update.cpp	2008-04-10 14:53:37.000000000 -0600
@@ -102,6 +102,7 @@
     force->nktv2p = 1.0;
     force->qqr2e = 1.0;
     force->qe2f = 1.0;
+    force->vxmu2f = 1.0;
     dt = 0.005;
     neighbor->skin = 0.3;
     
@@ -112,6 +113,7 @@
     force->nktv2p = 68568.415;
     force->qqr2e = 332.06371;
     force->qe2f = 23.060549; 
+    force->vxmu2f = 1.4393264316e4;
     dt = 1.0;
     neighbor->skin = 2.0;
 
@@ -122,6 +124,7 @@
     force->nktv2p = 1.6021765e6;
     force->qqr2e = 14.399645;
     force->qe2f = 1.0;
+    force->vxmu2f = 0.6241509647;
     dt = 0.001;
     neighbor->skin = 2.0;
   } else error->all("Illegal units command");