diff --git a/doc/Eqs/fix_nphug.jpg b/doc/Eqs/fix_nphug.jpg
new file mode 100644
index 000000000..beed5ed5c
Binary files /dev/null and b/doc/Eqs/fix_nphug.jpg differ
diff --git a/doc/Eqs/fix_nphug.tex b/doc/Eqs/fix_nphug.tex
new file mode 100644
index 000000000..4e162b69b
--- /dev/null
+++ b/doc/Eqs/fix_nphug.tex
@@ -0,0 +1,9 @@
+\documentstyle[12pt]{article}
+
+\begin{document}
+
+$$
+T_t - T = \frac{\left(\frac{1}{2}\left(P + P_0\right)\left(V_0 - V\right) + E_0 - E\right)}{N_{dof} k_B } = Delta
+$$
+
+\end{document}
diff --git a/doc/fix.txt b/doc/fix.txt
index 23a478c49..428b377fa 100644
--- a/doc/fix.txt
+++ b/doc/fix.txt
@@ -1,271 +1,272 @@
 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
 
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 
 :line
 
 fix command :h3
 
 [Syntax:]
 
 fix ID group-ID style args :pre
 
 ID = user-assigned name for the fix
 group-ID = ID of the group of atoms to apply the fix to
 style = one of a long list of possible style names (see below)
 args = arguments used by a particular style :ul
 
 [Examples:]
 
 fix 1 all nve
 fix 3 all nvt temp 300.0 300.0 0.01
 fix mine top setforce 0.0 NULL 0.0 :pre
 
 [Description:]
 
 Set a fix that will be applied to a group of atoms.  In LAMMPS, a
 "fix" is any operation that is applied to the system during
 timestepping or minimization.  Examples include updating of atom
 positions and velocities due to time integration, controlling
 temperature, applying constraint forces to atoms, enforcing boundary
 conditions, computing diagnostics, etc.  There are dozens of fixes
 defined in LAMMPS and new ones can be added; see "this
 section"_Section_modify.html for a discussion.
 
 Fixes perform their operations at different stages of the timestep.
 If 2 or more fixes operate at the same stage of the timestep, they are
 invoked in the order they were specified in the input script.
 
 The ID of a fix can only contain alphanumeric characters and
 underscores.
 
 Fixes can be deleted with the "unfix"_unfix.html command.
 
 IMPORTANT NOTE: The "unfix"_unfix.html command is the only way to turn
 off a fix; simply specifying a new fix with a similar style will not
 turn off the first one.  This is especially important to realize for
 integration fixes.  For example, using a "fix nve"_fix_nve.html
 command for a second run after using a "fix nvt"_fix_nh.html command
 for the first run, will not cancel out the NVT time integration
 invoked by the "fix nvt" command.  Thus two time integrators would be
 in place!
 
 If you specify a new fix with the same ID and style as an existing
 fix, the old fix is deleted and the new one is created (presumably
 with new settings).  This is the same as if an "unfix" command were
 first performed on the old fix, except that the new fix is kept in the
 same order relative to the existing fixes as the old one originally
 was.  Note that this operation also wipes out any additional changes
 made to the old fix via the "fix_modify"_fix_modify.html command.
 
 The "fix modify"_fix_modify.html command allows settings for some
 fixes to be reset.  See the doc page for individual fixes for details.
 
 Some fixes store an internal "state" which is written to binary
 restart files via the "restart"_restart.html or
 "write_restart"_write_restart.html commands.  This allows the fix to
 continue on with its calculations in a restarted simulation.  See the
 "read_restart"_read_restart.html command for info on how to re-specify
 a fix in an input script that reads a restart file.  See the doc pages
 for individual fixes for info on which ones can be restarted.
 
 :line
 
 Some fixes calculate one of three styles of quantities: global,
 per-atom, or local, which can be used by other commands or output as
 described below.  A global quantity is one or more system-wide values,
 e.g. the energy of a wall interacting with particles.  A per-atom
 quantity is one or more values per atom, e.g. the displacement vector
 for each atom since time 0.  Per-atom values are set to 0.0 for atoms
 not in the specified fix group.  Local quantities are calculated by
 each processor based on the atoms it owns, but there may be zero or
 more per atoms.
 
 Note that a single fix may produces either global or per-atom or local
 quantities (or none at all), but never more than one of these.
 
 Global, per-atom, and local quantities each come in three kinds: a
 single scalar value, a vector of values, or a 2d array of values.  The
 doc page for each fix describes the style and kind of values it
 produces, e.g. a per-atom vector.  Some fixes produce more than one
 kind of a single style, e.g. a global scalar and a global vector.
 
 When a fix quantity is accessed, as in many of the output commands
 discussed below, it can be referenced via the following bracket
 notation, where ID is the ID of the fix:
 
 f_ID | entire scalar, vector, or array
 f_ID\[I\] | one element of vector, one column of array
 f_ID\[I\]\[J\] | one element of array :tb(s=|)
 
 In other words, using one bracket reduces the dimension of the
 quantity once (vector -> scalar, array -> vector).  Using two brackets
 reduces the dimension twice (array -> scalar).  Thus a command that
 uses scalar fix values as input can also process elements of a vector
 or array.
 
 Note that commands and "variables"_variable.html which use fix
 quantities typically do not allow for all kinds, e.g. a command may
 require a vector of values, not a scalar.  This means there is no
 ambiguity about referring to a fix quantity as f_ID even if it
 produces, for example, both a scalar and vector.  The doc pages for
 various commands explain the details.
 
 :line
 
 In LAMMPS, the values generated by a fix can be used in several ways:
 
 Global values can be output via the "thermo_style
 custom"_thermo_style.html or "fix ave/time"_fix_ave_time.html command.
 Or the values can be referenced in a "variable equal"_variable.html or
 "variable atom"_variable.html command. :ulb,l
 
 Per-atom values can be output via the "dump custom"_dump.html command
 or the "fix ave/spatial"_fix_ave_spatial.html command.  Or they can be
 time-averaged via the "fix ave/atom"_fix_ave_atom.html command or
 reduced by the "compute reduce"_compute_reduce.html command.  Or the
 per-atom values can be referenced in an "atom-style
 variable"_variable.html. :l
 
 Local values can be reduced by the "compute
 reduce"_compute_reduce.html command, or histogrammed by the "fix
 ave/histo"_fix_ave_histo.html command. :l,ule
 
 See this "howto section"_Section_howto.html#howto_15 for a summary of
 various LAMMPS output options, many of which involve fixes.
 
 The results of fixes that calculate global quantities can be either
 "intensive" or "extensive" values.  Intensive means the value is
 independent of the number of atoms in the simulation,
 e.g. temperature.  Extensive means the value scales with the number of
 atoms in the simulation, e.g. total rotational kinetic energy.
 "Thermodynamic output"_thermo_style.html will normalize extensive
 values by the number of atoms in the system, depending on the
 "thermo_modify norm" setting.  It will not normalize intensive values.
 If a fix value is accessed in another way, e.g. by a
 "variable"_variable.html, you may want to know whether it is an
 intensive or extensive value.  See the doc page for individual fixes
 for further info.
 
 :line
 
 Each fix style has its own documentation page which describes its
 arguments and what it does, as listed below.  Here is an alphabetic
 list of fix styles available in LAMMPS:
 
 "adapt"_fix_adapt.html - change a simulation parameter over time
 "addforce"_fix_addforce.html - add a force to each atom
 "aveforce"_fix_aveforce.html - add an averaged force to each atom
 "ave/atom"_fix_ave_atom.html - compute per-atom time-averaged quantities
 "ave/histo"_fix_ave_atom.html - compute/output time-averaged histograms
 "ave/spatial"_fix_ave_spatial.html - compute/output time-averaged per-atom quantities by layer
 "ave/time"_fix_ave_time.html - compute/output global time-averaged quantities
 "bond/break"_fix_bond_break.html - break bonds on the fly
 "bond/create"_fix_bond_create.html - create bonds on the fly
 "bond/swap"_fix_bond_swap.html - Monte Carlo bond swapping
 "box/relax"_fix_box_relax.html - relax box size during energy minimization
 "deform"_fix_deform.html - change the simulation box size/shape
 "deposit"_fix_deposit.html - add new atoms above a surface
 "drag"_fix_drag.html - drag atoms towards a defined coordinate
 "dt/reset"_fix_dt_reset.html - reset the timestep based on velocity, forces
 "efield"_fix_efield.html - impose electric field on system
 "enforce2d"_fix_enforce2d.html - zero out z-dimension velocity and force
 "evaporate"_fix_evaporate.html - remove atoms from simulation periodically
 "external"_fix_external.html - callback to an external driver program
 "freeze"_fix_freeze.html - freeze atoms in a granular simulation
 "gravity"_fix_gravity.html - add gravity to atoms in a granular simulation
 "gcmc"_fix_gcmc.html - grand canonical insertions/deletions
 "heat"_fix_heat.html - add/subtract momentum-conserving heat
 "indent"_fix_indent.html - impose force due to an indenter
 "langevin"_fix_langevin.html - Langevin temperature control
 "lineforce"_fix_lineforce.html - constrain atoms to move in a line
 "momentum"_fix_momentum.html - zero the linear and/or angular momentum of a group of atoms
 "move"_fix_move.html - move atoms in a prescribed fashion
 "msst"_fix_msst.html - multi-scale shock technique (MSST) integration
 "neb"_fix_neb.html - nudged elastic band (NEB) spring forces
 "nph"_fix_nh.html - constant NPH time integration via Nose/Hoover
 "nph/asphere"_fix_nph_asphere.html - NPH for aspherical particles
 "nph/sphere"_fix_nph_sphere.html - NPH for spherical particles
+"nphug"_fix_nphug.html - Constant-stress Hugoniostat integration
 "npt"_fix_nh.html - constant NPT time integration via Nose/Hoover
 "npt/asphere"_fix_npt_asphere.html - NPT for aspherical particles
 "npt/sphere"_fix_npt_sphere.html - NPT for spherical particles
 "nve"_fix_nve.html - constant NVE time integration
 "nve/asphere"_fix_nve_asphere.html - NVT for aspherical particles
 "nve/limit"_fix_nve_limit.html - NVE with limited step length
 "nve/noforce"_fix_nve_noforce.html - NVE without forces (v only)
 "nve/sphere"_fix_nve_sphere.html - NVT for spherical particles
 "nvt"_fix_nh.html - constant NVT time integration via Nose/Hoover
 "nvt/asphere"_fix_nvt_asphere.html - NVT for aspherical particles
 "nvt/sllod"_fix_nvt_sllod.html - NVT for NEMD with SLLOD equations
 "nvt/sphere"_fix_nvt_sphere.html - NVT for spherical particles
 "orient/fcc"_fix_orient_fcc.html - add grain boundary migration force
 "planeforce"_fix_planeforce.html - constrain atoms to move in a plane
 "poems"_fix_poems.html - constrain clusters of atoms to move \
   as coupled rigid bodies
 "pour"_fix_pour.html - pour new atoms into a granular simulation domain
 "press/berendsen"_fix_press_berendsen.html - pressure control by \
      Berendsen barostat
 "print"_fix_print.html - print text and variables during a simulation
 "reax/bonds"_fix_reax_bonds.html - write out ReaxFF bond information \
 "recenter"_fix_recenter.html - constrain the center-of-mass position \
   of a group of atoms
 "restrain"_fix_restrain.html - constrain a bond, angle, dihedral
 "rigid"_fix_rigid.html - constrain one or more clusters of atoms to \
      move as a rigid body with NVE integration
 "rigid/nve"_fix_rigid.html - constrain one or more clusters of atoms to \
      move as a rigid body with alternate NVE integration
 "rigid/nvt"_fix_rigid.html - constrain one or more clusters of atoms to \
      move as a rigid body with NVT integration
 "setforce"_fix_setforce.html - set the force on each atom
 "shake"_fix_shake.html - SHAKE constraints on bonds and/or angles
 "spring"_fix_spring.html - apply harmonic spring force to group of atoms
 "spring/rg"_fix_spring_rg.html - spring on radius of gyration of \
      group of atoms
 "spring/self"_fix_spring_self.html - spring from each atom to its origin
 "srd"_fix_srd.html - stochastic rotation dynamics (SRD)
 "store/force"_fix_store_force.html - store force on each atom
 "store/state"_fix_store_state.html - store attributes for each atom
 "temp/berendsen"_fix_temp_berendsen.html - temperature control by \
      Berendsen thermostat
 "temp/rescale"_fix_temp_rescale.html - temperature control by \
      velocity rescaling
 "thermal/conductivity"_fix_thermal_conductivity.html - Muller-Plathe kinetic energy exchange for \
      thermal conductivity calculation
 "tmd"_fix_tmd.html - guide a group of atoms to a new configuration
 "ttm"_fix_ttm.html - two-temperature model for electronic/atomic coupling
 "viscosity"_fix_viscosity.html - Muller-Plathe momentum exchange for \
      viscosity calculation
 "viscous"_fix_viscous.html - viscous damping for granular simulations
 "wall/colloid"_fix_wall.html - Lennard-Jones wall interacting with finite-size particles
 "wall/gran"_fix_wall_gran.html - frictional wall(s) for granular simulations
 "wall/harmonic"_fix_wall.html - harmonic spring wall
 "wall/lj126"_fix_wall.html - Lennard-Jones 12-6 wall
 "wall/lj93"_fix_wall.html - Lennard-Jones 9-3 wall
 "wall/reflect"_fix_wall_reflect.html - reflecting wall(s)
 "wall/region"_fix_wall_region.html - use region surface as wall
 "wall/srd"_fix_wall_srd.html - slip/no-slip wall for SRD particles :ul
 
 There are also additional fix styles submitted by users which are
 included in the LAMMPS distribution.  The list of these with links to
 the individual styles are given in the fix section of "this
 page"_Section_commands.html#cmd_5.
 
 There are also additional accelerated fix styles included in the
 LAMMPS distribution for faster performance on CPUs and GPUs.  The list
 of these with links to the individual styles are given in the pair
 section of "this page"_Section_commands.html#cmd_5.
 
 [Restrictions:]
 
 Some fix styles are part of specific packages.  They are only enabled
 if LAMMPS was built with that package.  See the "Making
 LAMMPS"_Section_start.html#start_3 section for more info on packages.
 The doc pages for individual fixes tell if it is part of a package.
 
 [Related commands:]
 
 "unfix"_unfix.html, "fix_modify"_fix_modify.html
 
 [Default:] none
diff --git a/doc/fix_nphug.html b/doc/fix_nphug.html
index e4ad3a3f2..5d4b56cbc 100644
--- a/doc/fix_nphug.html
+++ b/doc/fix_nphug.html
@@ -1,219 +1,214 @@
 <HTML>
 <CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> 
 </CENTER>
 
 
 
 
 
 
 <HR>
 
 <H3>fix nphug command 
 </H3>
 <P><B>Syntax:</B>
 </P>
 <PRE>fix ID group-ID nphug keyword value ... 
 </PRE>
 <UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command 
 
 <PRE>one or more keyword value pairs may be appended
 keyword = <I>temp</I> or <I>iso</I> or <I>aniso</I> or <I>tri</I> or <I>x</I> or <I>y</I> or <I>z</I> or <I>couple</I> or <I>tchain</I> or <I>pchain</I> or <I>mtk</I> or <I>tloop</I> or <I>ploop</I> or <I>nreset</I> or <I>drag</I> or <I>dilate</I> or <I>scaleyz</I> or <I>scalexz</I> or <I>scalexy</I>
   <I>temp</I> values = Value1 Value2 Tdamp
     Value1, Value2 = Nose-Hoover target temperatures, ignored by Hugoniostat
     Tdamp = temperature damping parameter (time units)
   <I>iso</I> or <I>aniso</I> or <I>tri</I> values = Pstart Pstop Pdamp
     Pstart,Pstop = scalar external pressures, must be equal (pressure units)
     Pdamp = pressure damping parameter (time units)
   <I>x</I> or <I>y</I> or <I>z</I> or <I>xy</I> or <I>yz</I> or <I>xz</I> values = Pstart Pstop Pdamp
     Pstart,Pstop = external stress tensor components, must be equal (pressure units)
     Pdamp = stress damping parameter (time units)
   <I>couple</I> = <I>none</I> or <I>xyz</I> or <I>xy</I> or <I>yz</I> or <I>xz</I>
   <I>tchain</I> value = length of thermostat chain (1 = single thermostat)
   <I>pchain</I> values = length of thermostat chain on barostat (0 = no thermostat)
   <I>mtk</I> value = <I>yes</I> or <I>no</I> = add in MTK adjustment term or not
   <I>tloop</I> value = number of sub-cycles to perform on thermostat
   <I>ploop</I> value = number of sub-cycles to perform on barostat thermostat
   <I>nreset</I> value = reset reference cell every this many timesteps
   <I>drag</I> value = drag factor added to barostat/thermostat (0.0 = no drag)
   <I>dilate</I> value = <I>all</I> or <I>partial</I>
   <I>scaleyz</I> value = <I>yes</I> or <I>no</I> = scale yz with lz
   <I>scalexz</I> value = <I>yes</I> or <I>no</I> = scale xz with lz
   <I>scalexy</I> value = <I>yes</I> or <I>no</I> = scale xy with ly 
 </PRE>
 
 </UL>
 <P><B>Examples:</B>
 </P>
-<P>fix myhug all nphug temp 1.0 1.0 10.0 z 40.0 40.0 70.0
+<PRE>fix myhug all nphug temp 1.0 1.0 10.0 z 40.0 40.0 70.0
 fix myhug all nphug temp 1.0 1.0 10.0 iso 40.0 40.0 70.0 drag 200.0 tchain 1 pchain 0 
-</P>
+</PRE>
 <P><B>Description:</B>
 </P>
 <P>This command is a variant of the Nose-Hoover
 <A HREF = "fix_nh.html">fix npt</A> fix style.
 It performs time integration of the Hugoniostat equations
 of motion developed by Ravelo et al. <A HREF = "#Ravelo">(Ravelo)</A>.
 These equations compress the system to a state with average
 axial stress or pressure equal to the specified target value 
 and that satisfies the Rankine-Hugoniot (RH)
 jump conditions for steady shocks.
 </P>
 <P>The compression can be performed 
 either
-hydrostatically (using keyword iso, aniso, or tri) or uniaxially
-(using keywords x, y, or z).  In the hydrostatic case,
+hydrostatically (using keyword <I>iso</I>, <I>aniso</I>, or <I>tri</I>) or uniaxially
+(using keywords <I>x</I>, <I>y</I>, or <I>z</I>).  In the hydrostatic case,
 the cell dimensions change dynamically so that the average axial stress
 in all three directions converges towards the specified target value. 
 In the uniaxial case, the chosen cell dimension changes dynamically 
 so that the average
 axial stress in that direction converges towards the target value. The
 other two cell dimensions are kept fixed (zero lateral strain).
 </P>
 <P>This leads to the following additional restrictions on the keywords:
 </P>
-<P>One and only one of the following keywords should be used: iso, aniso, tri, x, y, z, 
-</P>
-<P>The specified initial and final target pressures must be the same.
-</P>
-<P>The keywords xy, xz, yz may not be used.
-</P>
-<P>The only admissible value for the couple keyword is xyz, 
-which has the same effect as keyword iso 
-</P>
-<UL><LI>The temp keyword must be used in order to specify the time constant for 
-<LI>kinetic energy relaxation, but initial and final target temperature values
-<LI>are ignored.  
+<UL><LI>One and only one of the following keywords should be used: <I>iso</I>, <I>aniso</I>, <I>tri</I>, <I>x</I>, <I>y</I>, <I>z</I> 
+<LI>The specified initial and final target pressures must be the same.
+<LI>The keywords <I>xy</I>, <I>xz</I>, <I>yz</I> may not be used.
+<LI>The only admissible value for the couple keyword is <I>xyz</I>, which has the same effect as keyword <I>iso</I> 
+<LI>The <I>temp</I> keyword must be used to specify the time constant for kinetic energy relaxation, but initial and final target temperature values are ignored.  
 </UL>
 <P>Essentially, a Hugoniostat simulation is an NPT simulation in which the
 user-specified target temperature is replaced with a time-dependent 
 target temperature Tt obtained from the following equation:
 </P>
-<P>Tt - T = (0.5*(P+P0)(V0-V)+E0-E)/(Ndof kB) = Delta
-</P>
+<CENTER><IMG SRC = "Eqs/fix_nphug.jpg">
+</CENTER>
 <P>where T and Tt are the instantaneous and target temperatures,
-P and P0 are the instantaneous and reference pressure or axial stress,
+P and P0 are the instantaneous and reference pressures or axial stresses,
 depending on whether hydrostatic or uniaxial compression is being 
 performed, V and V0 are the instantaneous and reference volumes,
 E and E0 are the instantaneous and reference internal energy (potential
 plus kinetic), Ndof is the number of degrees of freedom used in the
 definition of temperature, and kB is the Boltzmann constant. Delta is the 
 negative deviation of the instantaneous temperature from the target temperature.
-The values of E0, V0, and P0 are the instantaneous values at the start of
+When the system reaches a stable equilibrium, the value of Delta should 
+fluctuate about zero.
+</P>
+<P>The values of E0, V0, and P0 are the instantaneous values at the start of
 the simulation. These can be overridden using the fix_modify keywords <I>e0</I>,
 <I>v0</I>, and <I>p0</I> described below.
 </P>
 <HR>
 
 <P>IMPORTANT NOTE: Unlike the <A HREF = "fix_temp_berendsen.html">fix
 temp/berendsen</A> command which performs
 thermostatting but NO time integration, this fix performs
 thermostatting/barostatting AND time integration.  Thus you should not
 use any other time integration fix, such as <A HREF = "fix_nve.html">fix nve</A> on
 atoms to which this fix is applied.  Likewise, this fix should not be 
 used on atoms that have their temperature controlled by another fix 
 - e.g. by <A HREF = "fix_nh.html">fix langevin</A> or <A HREF = "fix_temp_rescale.html">fix temp/rescale</A>
 commands.
 </P>
 <HR>
 
-<P>This fix compute a temperature and pressure each timestep.  To do
+<P>This fix computes a temperature and pressure at each timestep.  To do
 this, the fix creates its own computes of style "temp" and "pressure",
 as if one of these two sets of commands had been issued:
 </P>
 <PRE>compute fix-ID_temp group-ID temp
 compute fix-ID_press group-ID pressure fix-ID_temp 
 </PRE>
 <PRE>compute fix-ID_temp all temp
 compute fix-ID_press all pressure fix-ID_temp 
 </PRE>
 <P>See the <A HREF = "compute_temp.html">compute temp</A> and <A HREF = "compute_pressure.html">compute
 pressure</A> commands for details.  Note that the
 IDs of the new computes are the fix-ID + underscore + "temp" or fix_ID
-+ underscore + "press".  For fix nvt, the group for the new computes
-is the same as the fix group.  For fix nph and fix npt, the group for
++ underscore + "press".  The group for
 the new computes is "all" since pressure is computed for the entire
 system.
 </P>
 <P>Note that these are NOT the computes used by thermodynamic output (see
 the <A HREF = "thermo_style.html">thermo_style</A> command) with ID = <I>thermo_temp</I>
 and <I>thermo_press</I>.  This means you can change the attributes of this
 fix's temperature or pressure via the
 <A HREF = "compute_modify.html">compute_modify</A> command or print this temperature
 or pressure during thermodynamic output via the <A HREF = "thermo_style.html">thermo_style
 custom</A> command using the appropriate compute-ID.
 It also means that changing attributes of <I>thermo_temp</I> or
 <I>thermo_press</I> will have no effect on this fix.
 </P>
 <HR>
 
 <P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
 </P>
 <P>This fix writes the values of E0, V0, and P0, as well as the 
 state of all the thermostat and barostat
 variables to <A HREF = "restart.html">binary restart files</A>.  See the
 <A HREF = "read_restart.html">read_restart</A> command for info on how to re-specify
 a fix in an input script that reads a restart file, so that the
 operation of the fix continues in an uninterrupted fashion.
 </P>
 <P>The <A HREF = "fix_modify.html">fix_modify</A> <I>e0</I>, <I>v0</I> and <I>p0</I> keywords
-can be used to define the values of these quantities. Note the
-the values for <I>e0</I> and <I>v0</I> are extensive, and so must correspond
+can be used to define the values of E0, V0, and P0. Note the
+the values for <I>e0</I> and <I>p0</I> are extensive, and so must correspond
 to the total energy and volume of the entire system, not energy and 
 volume per atom. If any of these quantities are not specified, then the
 instanteous value in the system at the start of the simulation is used.
 </P>
 <P>The <A HREF = "fix_modify.html">fix_modify</A> <I>temp</I> and <I>press</I> options are
 supported by these fixes.  You can use them to assign a
 <A HREF = "compute.html">compute</A> you have defined to this fix which will be used
 in its thermostatting or barostatting procedure, as described above.
 If you do this, note that the kinetic energy derived from the compute
 temperature should be consistent with the virial term computed using
 all atoms for the pressure.  LAMMPS will warn you if you choose to
 compute temperature on a subset of atoms.
 </P>
 <P>The <A HREF = "fix_modify.html">fix_modify</A> <I>energy</I> option is supported by these
 fixes to add the energy change induced by Nose/Hoover thermostatting
 and barostatting to the system's potential energy as part of
 <A HREF = "thermo_style.html">thermodynamic output</A>. Either way, this energy is *not*
 included in the definition of internal energy E when calculating the value
 of Delta in the above equation.
 </P>
 <P>These fixes compute a global scalar and a global vector of quantities,
 which can be accessed by various <A HREF = "Section_howto.html#howto_15">output
 commands</A>.  The scalar value calculated by
 these fixes is "extensive"; the vector values are "intensive".
 </P>
 <P>The scalar is the cumulative energy change due to the fix.
 </P>
 <P>The vector stores three quantities unique to this fix (Delta, Us, and up),
 followed by all the internal Nose/Hoover thermostat and barostat
 variables defined for <A HREF = "fix_nh.html">fix_style npt</A>. Delta is the deviation
 of the temperature from the target temperature, given by the above equation.
-Us and up are the shock and particle velocity correspponding to a steady
-shock calculated from the Rh conditions. They have units of distance/time.
+Us and up are the shock and particle velocity corresponding to a steady
+shock calculated from the RH conditions. They have units of distance/time.
 </P>
 <P><B>Restrictions:</B>
 </P>
 <P>This fix style is part of the SHOCK package.  It is only enabled if
 LAMMPS was built with that package. See the <A HREF = "Section_start.html#start_3">Making
 LAMMPS</A> section for more info.
 </P>
 <P>All the usual restrictions for <A HREF = "fix_nh.html">fix_style npt</A> apply,
 plus the additional ones mentioned above.
 </P>
 <P><B>Related commands:</B>
 </P>
 <P><A HREF = "fix_msst.html">fix msst</A>, <A HREF = "fix_nh.html">fix npt</A>, <A HREF = "fix_modify.html">fix_modify</A>
 </P>
 <P><B>Default:</B>
 </P>
 <P>The keyword defaults are the same as those for <A HREF = "fix_nh.html">fix npt</A>
 </P>
 <HR>
 
 <A NAME = "Ravelo"></A>
 
 <P><B>(Ravelo)</B> Ravelo, Holian, Germann and Lomdahl, Phys Rev B, 70, 014103 (2004).
 </P>
 </HTML>
diff --git a/doc/fix_nphug.txt b/doc/fix_nphug.txt
index 6bb566430..6ac59d1e1 100644
--- a/doc/fix_nphug.txt
+++ b/doc/fix_nphug.txt
@@ -1,211 +1,206 @@
 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
 
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 
 :line
 
 fix nphug command :h3
 
 [Syntax:]
 
 fix ID group-ID nphug keyword value ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 one or more keyword value pairs may be appended
 keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scaleyz} or {scalexz} or {scalexy}
   {temp} values = Value1 Value2 Tdamp
     Value1, Value2 = Nose-Hoover target temperatures, ignored by Hugoniostat
     Tdamp = temperature damping parameter (time units)
   {iso} or {aniso} or {tri} values = Pstart Pstop Pdamp
     Pstart,Pstop = scalar external pressures, must be equal (pressure units)
     Pdamp = pressure damping parameter (time units)
   {x} or {y} or {z} or {xy} or {yz} or {xz} values = Pstart Pstop Pdamp
     Pstart,Pstop = external stress tensor components, must be equal (pressure units)
     Pdamp = stress damping parameter (time units)
   {couple} = {none} or {xyz} or {xy} or {yz} or {xz}
   {tchain} value = length of thermostat chain (1 = single thermostat)
   {pchain} values = length of thermostat chain on barostat (0 = no thermostat)
   {mtk} value = {yes} or {no} = add in MTK adjustment term or not
   {tloop} value = number of sub-cycles to perform on thermostat
   {ploop} value = number of sub-cycles to perform on barostat thermostat
   {nreset} value = reset reference cell every this many timesteps
   {drag} value = drag factor added to barostat/thermostat (0.0 = no drag)
   {dilate} value = {all} or {partial}
   {scaleyz} value = {yes} or {no} = scale yz with lz
   {scalexz} value = {yes} or {no} = scale xz with lz
   {scalexy} value = {yes} or {no} = scale xy with ly :pre
 :ule
 
 [Examples:]
 
 fix myhug all nphug temp 1.0 1.0 10.0 z 40.0 40.0 70.0
-fix myhug all nphug temp 1.0 1.0 10.0 iso 40.0 40.0 70.0 drag 200.0 tchain 1 pchain 0 
+fix myhug all nphug temp 1.0 1.0 10.0 iso 40.0 40.0 70.0 drag 200.0 tchain 1 pchain 0 :pre
 
 [Description:]
 
 This command is a variant of the Nose-Hoover
 "fix npt"_fix_nh.html fix style.
 It performs time integration of the Hugoniostat equations
 of motion developed by Ravelo et al. "(Ravelo)"_#Ravelo.
 These equations compress the system to a state with average
 axial stress or pressure equal to the specified target value 
 and that satisfies the Rankine-Hugoniot (RH)
 jump conditions for steady shocks.
 
 The compression can be performed 
 either
-hydrostatically (using keyword iso, aniso, or tri) or uniaxially
-(using keywords x, y, or z).  In the hydrostatic case,
+hydrostatically (using keyword {iso}, {aniso}, or {tri}) or uniaxially
+(using keywords {x}, {y}, or {z}).  In the hydrostatic case,
 the cell dimensions change dynamically so that the average axial stress
 in all three directions converges towards the specified target value. 
 In the uniaxial case, the chosen cell dimension changes dynamically 
 so that the average
 axial stress in that direction converges towards the target value. The
 other two cell dimensions are kept fixed (zero lateral strain).
 
 This leads to the following additional restrictions on the keywords:
 
-One and only one of the following keywords should be used: iso, aniso, tri, x, y, z, 
-
+One and only one of the following keywords should be used: {iso}, {aniso}, {tri}, {x}, {y}, {z} 
 The specified initial and final target pressures must be the same.
-
-The keywords xy, xz, yz may not be used.
-
-The only admissible value for the couple keyword is xyz, 
-which has the same effect as keyword iso 
-
-The temp keyword must be used in order to specify the time constant for 
-kinetic energy relaxation, but initial and final target temperature values
-are ignored.  :ul
+The keywords {xy}, {xz}, {yz} may not be used.
+The only admissible value for the couple keyword is {xyz}, which has the same effect as keyword {iso} 
+The {temp} keyword must be used to specify the time constant for kinetic energy relaxation, but initial and final target temperature values are ignored.  :ul
 
 Essentially, a Hugoniostat simulation is an NPT simulation in which the
 user-specified target temperature is replaced with a time-dependent 
 target temperature Tt obtained from the following equation:
 
-Tt - T = (0.5*(P+P0)(V0-V)+E0-E)/(Ndof kB) = Delta
+:c,image(Eqs/fix_nphug.jpg) 
 
 where T and Tt are the instantaneous and target temperatures,
-P and P0 are the instantaneous and reference pressure or axial stress,
+P and P0 are the instantaneous and reference pressures or axial stresses,
 depending on whether hydrostatic or uniaxial compression is being 
 performed, V and V0 are the instantaneous and reference volumes,
 E and E0 are the instantaneous and reference internal energy (potential
 plus kinetic), Ndof is the number of degrees of freedom used in the
 definition of temperature, and kB is the Boltzmann constant. Delta is the 
 negative deviation of the instantaneous temperature from the target temperature.
+When the system reaches a stable equilibrium, the value of Delta should 
+fluctuate about zero.
+
 The values of E0, V0, and P0 are the instantaneous values at the start of
 the simulation. These can be overridden using the fix_modify keywords {e0},
 {v0}, and {p0} described below.
 
 :line
 
 IMPORTANT NOTE: Unlike the "fix
 temp/berendsen"_fix_temp_berendsen.html command which performs
 thermostatting but NO time integration, this fix performs
 thermostatting/barostatting AND time integration.  Thus you should not
 use any other time integration fix, such as "fix nve"_fix_nve.html on
 atoms to which this fix is applied.  Likewise, this fix should not be 
 used on atoms that have their temperature controlled by another fix 
 - e.g. by "fix langevin"_fix_nh.html or "fix temp/rescale"_fix_temp_rescale.html
 commands.
 
 :line
 
-This fix compute a temperature and pressure each timestep.  To do
+This fix computes a temperature and pressure at each timestep.  To do
 this, the fix creates its own computes of style "temp" and "pressure",
 as if one of these two sets of commands had been issued:
 
 compute fix-ID_temp group-ID temp
 compute fix-ID_press group-ID pressure fix-ID_temp :pre
 
 compute fix-ID_temp all temp
 compute fix-ID_press all pressure fix-ID_temp :pre
 
 See the "compute temp"_compute_temp.html and "compute
 pressure"_compute_pressure.html commands for details.  Note that the
 IDs of the new computes are the fix-ID + underscore + "temp" or fix_ID
-+ underscore + "press".  For fix nvt, the group for the new computes
-is the same as the fix group.  For fix nph and fix npt, the group for
++ underscore + "press".  The group for
 the new computes is "all" since pressure is computed for the entire
 system.
 
 Note that these are NOT the computes used by thermodynamic output (see
 the "thermo_style"_thermo_style.html command) with ID = {thermo_temp}
 and {thermo_press}.  This means you can change the attributes of this
 fix's temperature or pressure via the
 "compute_modify"_compute_modify.html command or print this temperature
 or pressure during thermodynamic output via the "thermo_style
 custom"_thermo_style.html command using the appropriate compute-ID.
 It also means that changing attributes of {thermo_temp} or
 {thermo_press} will have no effect on this fix.
 
 :line
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
 This fix writes the values of E0, V0, and P0, as well as the 
 state of all the thermostat and barostat
 variables to "binary restart files"_restart.html.  See the
 "read_restart"_read_restart.html command for info on how to re-specify
 a fix in an input script that reads a restart file, so that the
 operation of the fix continues in an uninterrupted fashion.
 
 The "fix_modify"_fix_modify.html {e0}, {v0} and {p0} keywords
-can be used to define the values of these quantities. Note the
-the values for {e0} and {v0} are extensive, and so must correspond
+can be used to define the values of E0, V0, and P0. Note the
+the values for {e0} and {p0} are extensive, and so must correspond
 to the total energy and volume of the entire system, not energy and 
 volume per atom. If any of these quantities are not specified, then the
 instanteous value in the system at the start of the simulation is used.
 
 The "fix_modify"_fix_modify.html {temp} and {press} options are
 supported by these fixes.  You can use them to assign a
 "compute"_compute.html you have defined to this fix which will be used
 in its thermostatting or barostatting procedure, as described above.
 If you do this, note that the kinetic energy derived from the compute
 temperature should be consistent with the virial term computed using
 all atoms for the pressure.  LAMMPS will warn you if you choose to
 compute temperature on a subset of atoms.
 
 The "fix_modify"_fix_modify.html {energy} option is supported by these
 fixes to add the energy change induced by Nose/Hoover thermostatting
 and barostatting to the system's potential energy as part of
 "thermodynamic output"_thermo_style.html. Either way, this energy is *not*
 included in the definition of internal energy E when calculating the value
 of Delta in the above equation.
 
 These fixes compute a global scalar and a global vector of quantities,
 which can be accessed by various "output
 commands"_Section_howto.html#howto_15.  The scalar value calculated by
 these fixes is "extensive"; the vector values are "intensive".
 
 The scalar is the cumulative energy change due to the fix.
 
 The vector stores three quantities unique to this fix (Delta, Us, and up),
 followed by all the internal Nose/Hoover thermostat and barostat
 variables defined for "fix_style npt"_fix_nh.html. Delta is the deviation
 of the temperature from the target temperature, given by the above equation.
-Us and up are the shock and particle velocity correspponding to a steady
-shock calculated from the Rh conditions. They have units of distance/time.
+Us and up are the shock and particle velocity corresponding to a steady
+shock calculated from the RH conditions. They have units of distance/time.
 
 [Restrictions:]
 
 This fix style is part of the SHOCK package.  It is only enabled if
 LAMMPS was built with that package. See the "Making
 LAMMPS"_Section_start.html#start_3 section for more info.
 
 All the usual restrictions for "fix_style npt"_fix_nh.html apply,
 plus the additional ones mentioned above.
 
 [Related commands:]
 
 "fix msst"_fix_msst.html, "fix npt"_fix_nh.html, "fix_modify"_fix_modify.html
 
 [Default:]
 
 The keyword defaults are the same as those for "fix npt"_fix_nh.html
 
 :line
 
 :link(Ravelo)
 [(Ravelo)] Ravelo, Holian, Germann and Lomdahl, Phys Rev B, 70, 014103 (2004).