diff --git a/doc/Section_commands.html b/doc/Section_commands.html
index e758583e7..80faeefde 100644
--- a/doc/Section_commands.html
+++ b/doc/Section_commands.html
@@ -1,484 +1,484 @@
 <HTML>
 <CENTER><A HREF = "Section_start.html">Previous Section</A> - <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> - <A HREF = "Section_howto.html">Next Section</A> 
 </CENTER>
 
 
 
 
 
 
 <HR>
 
 <H3>3. Commands 
 </H3>
 <P>This section describes how a LAMMPS input script is formatted and what
 commands are used to define a LAMMPS simulation.
 </P>
 3.1 <A HREF = "#3_1">LAMMPS input script</A><BR>
 3.2 <A HREF = "#3_2">Parsing rules</A><BR>
 3.3 <A HREF = "#3_3">Input script structure</A><BR>
 3.4 <A HREF = "#3_4">Commands listed by category</A><BR>
 3.5 <A HREF = "#3_5">Commands listed alphabetically</A> <BR>
 
 <HR>
 
 <A NAME = "3_1"></A><H4>3.1 LAMMPS input script 
 </H4>
 <P>LAMMPS executes by reading commands from a input script (text file),
 one line at a time.  When the input script ends, LAMMPS exits.  Each
 command causes LAMMPS to take some action.  It may set an internal
 variable, read in a file, or run a simulation.  Most commands have
 default settings, which means you only need to use the command if you
 wish to change the default.
 </P>
 <P>In many cases, the ordering of commands in an input script is not
 important.  However the following rules apply:
 </P>
 <P>(1) LAMMPS does not read your entire input script and then perform a
 simulation with all the settings.  Rather, the input script is read
 one line at a time and each command takes effect when it is read.
 Thus this sequence of commands:
 </P>
 <PRE>timestep 0.5 
 run      100 
 run      100 
 </PRE>
 <P>does something different than this sequence:
 </P>
 <PRE>run      100 
 timestep 0.5 
 run      100 
 </PRE>
 <P>In the first case, the specified timestep (0.5 fmsec) is used for two
 simulations of 100 timesteps each.  In the 2nd case, the default
 timestep (1.0 fmsec) is used for the 1st 100 step simulation and a 0.5
 fmsec timestep is used for the 2nd one.
 </P>
 <P>(2) Some commands are only valid when they follow other commands.  For
 example you cannot set the temperature of a group of atoms until atoms
 have been defined and a group command is used to define which atoms
 belong to the group.
 </P>
 <P>(3) Sometimes command B will use values that can be set by command A.
 This means command A must precede command B in the input script if it
 is to have the desired effect.  For example, the
 <A HREF = "read_data.html">read_data</A> command initializes the system by setting
 up the simulation box and assigning atoms to processors.  If default
 values are not desired, the <A HREF = "processors.html">processors</A> and
 <A HREF = "boundary.html">boundary</A> commands need to be used before read_data to
 tell LAMMPS how to map processors to the simulation box.
 </P>
 <P>Many input script errors are detected by LAMMPS and an ERROR or
 WARNING message is printed.  <A HREF = "Section_errors.html">This section</A> gives
 more information on what errors mean.  The documentation for each
 command lists restrictions on how the command can be used.
 </P>
 <HR>
 
 <A NAME = "3_2"></A><H4>3.2 Parsing rules 
 </H4>
 <P>Each non-blank line in the input script is treated as a command.
 LAMMPS commands are case sensitive.  Command names are lower-case, as
 are specified command arguments.  Upper case letters may be used in
 file names or user-chosen ID strings.
 </P>
 <P>Here is how each line in the input script is parsed by LAMMPS:
 </P>
 <P>(1) If the last printable character on the line is a "&" character
 (with no surrounding quotes), the command is assumed to continue on
 the next line.  The next line is concatenated to the previous line by
 removing the "&" character and newline.  This allows long commands to
 be continued across two or more lines.
 </P>
 <P>(2) All characters from the first "#" character onward are treated as
 comment and discarded.  See an exception in (6).  Note that a
 comment after a trailing "&" character will prevent the command from
 continuing on the next line.  Also note that for multi-line commands a
 single leading "#" will comment out the entire command.
 </P>
 <P>(3) The line is searched repeatedly for $ characters, which indicate
 variables that are replaced with a text string.  See an exception in
 (6).  If the $ is followed by curly brackets, then the variable name
 is the text inside the curly brackets.  If no curly brackets follow
 the $, then the variable name is the single character immediately
 following the $.  Thus ${myTemp} and $x refer to variable names
 "myTemp" and "x".  See the <A HREF = "variable.html">variable</A> command for
 details of how strings are assigned to variables and how they are
 substituted for in input script commands.
 </P>
 <P>(4) The line is broken into "words" separated by whitespace (tabs,
 spaces).  Note that words can thus contain letters, digits,
 underscores, or punctuation characters.
 </P>
 <P>(5) The first word is the command name.  All successive words in the
 line are arguments.
 </P>
 <P>(6) If you want text with spaces to be treated as a single argument,
 it can be enclosed in double quotes.  The quotes are removed when the
 single argument is stored in its final form.  See the <A HREF = "dump_modify.html">dump modify
 format</A> or <A HREF = "if.html">if</A> commands for examples.  A "#"
 or "$" character that is between double quotes will not be treated as
 a comment indicator in (2) or substituted for as a variable in (3).
 </P>
 <HR>
 
 <H4><A NAME = "3_3"></A>3.3 Input script structure 
 </H4>
 <P>This section describes the structure of a typical LAMMPS input script.
 The "examples" directory in the LAMMPS distribution contains many
 sample input scripts; the corresponding problems are discussed in
 <A HREF = "Section_example.html">this section</A>, and animated on the <A HREF = "http://lammps.sandia.gov">LAMMPS WWW
 Site</A>.
 </P>
 <P>A LAMMPS input script typically has 4 parts:
 </P>
 <OL><LI>Initialization
 <LI>Atom definition
 <LI>Settings
 <LI>Run a simulation 
 </OL>
 <P>The last 2 parts can be repeated as many times as desired.  I.e. run a
 simulation, change some settings, run some more, etc.  Each of the 4
 parts is now described in more detail.  Remember that almost all the
 commands need only be used if a non-default value is desired.
 </P>
 <P>(1) Initialization
 </P>
 <P>Set parameters that need to be defined before atoms are created or
 read-in from a file.
 </P>
 <P>The relevant commands are <A HREF = "units.html">units</A>,
 <A HREF = "dimension.html">dimension</A>, <A HREF = "newton.html">newton</A>,
 <A HREF = "processors.html">processors</A>, <A HREF = "boundary.html">boundary</A>,
 <A HREF = "atom_style.html">atom_style</A>, <A HREF = "atom_modify.html">atom_modify</A>.
 </P>
 <P>If force-field parameters appear in the files that will be read, these
 commands tell LAMMPS what kinds of force fields are being used:
 <A HREF = "pair_style.html">pair_style</A>, <A HREF = "bond_style.html">bond_style</A>,
 <A HREF = "angle_style.html">angle_style</A>, <A HREF = "dihedral_style.html">dihedral_style</A>,
 <A HREF = "improper_style.html">improper_style</A>.
 </P>
 <P>(2) Atom definition
 </P>
 <P>There are 3 ways to define atoms in LAMMPS.  Read them in from a data
 or restart file via the <A HREF = "read_data.html">read_data</A> or
 <A HREF = "read_restart.html">read_restart</A> commands.  These files can contain
 molecular topology information.  Or create atoms on a lattice (with no
 molecular topology), using these commands: <A HREF = "lattice.html">lattice</A>,
 <A HREF = "region.html">region</A>, <A HREF = "create_box.html">create_box</A>,
 <A HREF = "create_atoms.html">create_atoms</A>.  The entire set of atoms can be
 duplicated to make a larger simulation using the
 <A HREF = "replicate.html">replicate</A> command.
 </P>
 <P>(3) Settings
 </P>
 <P>Once atoms and molecular topology are defined, a variety of settings
 can be specified: force field coefficients, simulation parameters,
 output options, etc.
 </P>
 <P>Force field coefficients are set by these commands (they can also be
 set in the read-in files): <A HREF = "pair_coeff.html">pair_coeff</A>,
 <A HREF = "bond_coeff.html">bond_coeff</A>, <A HREF = "angle_coeff.html">angle_coeff</A>,
 <A HREF = "dihedral_coeff.html">dihedral_coeff</A>,
 <A HREF = "improper_coeff.html">improper_coeff</A>,
 <A HREF = "kspace_style.html">kspace_style</A>, <A HREF = "dielectric.html">dielectric</A>,
 <A HREF = "special_bonds.html">special_bonds</A>.
 </P>
 <P>Various simulation parameters are set by these commands:
 <A HREF = "neighbor.html">neighbor</A>, <A HREF = "neigh_modify.html">neigh_modify</A>,
 <A HREF = "group.html">group</A>, <A HREF = "timestep.html">timestep</A>,
 <A HREF = "reset_timestep.html">reset_timestep</A>, <A HREF = "run_style.html">run_style</A>,
 <A HREF = "min_style.html">min_style</A>, <A HREF = "min_modify.html">min_modify</A>.
 </P>
 <P>Fixes impose a variety of boundary conditions, time integration, and
 diagnostic options.  The <A HREF = "fix.html">fix</A> command comes in many flavors.
 </P>
 <P>Various computations can be specified for execution during a
 simulation using the <A HREF = "compute.html">compute</A>,
 <A HREF = "compute_modify.html">compute_modify</A>, and <A HREF = "variable.html">variable</A>
 commands.
 </P>
 <P>Output options are set by the <A HREF = "thermo.html">thermo</A>, <A HREF = "dump.html">dump</A>,
 and <A HREF = "restart.html">restart</A> commands.
 </P>
 <P>(4) Run a simulation
 </P>
 <P>A molecular dynamics simulation is run using the <A HREF = "run.html">run</A>
 command.  Energy minimization (molecular statics) is performed using
 the <A HREF = "minimize.html">minimize</A> command.  A parallel tempering
 (replica-exchange) simulation can be run using the
 <A HREF = "temper.html">temper</A> command.
 </P>
 <HR>
 
 <A NAME = "3_4"></A><H4>3.4 Commands listed by category 
 </H4>
 <P>This section lists all LAMMPS commands, grouped by category.  The
 <A HREF = "#3_5">next section</A> lists the same commands alphabetically.  Note that
 some style options for some commands are part of specific LAMMPS
 packages, which means they cannot be used unless the package was
 included when LAMMPS was built.  Not all packages are included in a
 default LAMMPS build.  These dependencies are listed as Restrictions
 in the command's documentation.
 </P>
 <P>Initialization:
 </P>
 <P><A HREF = "atom_modify.html">atom_modify</A>, <A HREF = "atom_style.html">atom_style</A>,
 <A HREF = "boundary.html">boundary</A>, <A HREF = "dimension.html">dimension</A>,
 <A HREF = "newton.html">newton</A>, <A HREF = "processors.html">processors</A>, <A HREF = "units.html">units</A>
 </P>
 <P>Atom definition:
 </P>
 <P><A HREF = "create_atoms.html">create_atoms</A>, <A HREF = "create_box.html">create_box</A>,
 <A HREF = "lattice.html">lattice</A>, <A HREF = "read_data.html">read_data</A>,
 <A HREF = "read_restart.html">read_restart</A>, <A HREF = "region.html">region</A>,
 <A HREF = "replicate.html">replicate</A>
 </P>
 <P>Force fields:
 </P>
 <P><A HREF = "angle_coeff.html">angle_coeff</A>, <A HREF = "angle_style.html">angle_style</A>,
 <A HREF = "bond_coeff.html">bond_coeff</A>, <A HREF = "bond_style.html">bond_style</A>,
 <A HREF = "dielectric.html">dielectric</A>, <A HREF = "dihedral_coeff.html">dihedral_coeff</A>,
 <A HREF = "dihedral_style.html">dihedral_style</A>,
 <A HREF = "improper_coeff.html">improper_coeff</A>,
 <A HREF = "improper_style.html">improper_style</A>,
 <A HREF = "kspace_modify.html">kspace_modify</A>, <A HREF = "kspace_style.html">kspace_style</A>,
 <A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "pair_modify.html">pair_modify</A>,
 <A HREF = "pair_style.html">pair_style</A>, <A HREF = "pair_write.html">pair_write</A>,
 <A HREF = "special_bonds.html">special_bonds</A>
 </P>
 <P>Settings:
 </P>
 <P><A HREF = "communicate.html">communicate</A>, <A HREF = "dipole.html">dipole</A>,
 <A HREF = "group.html">group</A>, <A HREF = "mass.html">mass</A>, <A HREF = "min_modify.html">min_modify</A>,
 <A HREF = "min_style.html">min_style</A>, <A HREF = "neigh_modify.html">neigh_modify</A>,
 <A HREF = "neighbor.html">neighbor</A>, <A HREF = "reset_timestep.html">reset_timestep</A>,
 <A HREF = "run_style.html">run_style</A>, <A HREF = "set.html">set</A>, <A HREF = "shape.html">shape</A>,
 <A HREF = "timestep.html">timestep</A>, <A HREF = "velocity.html">velocity</A>
 </P>
 <P>Fixes:
 </P>
 <P><A HREF = "fix.html">fix</A>, <A HREF = "fix_modify.html">fix_modify</A>, <A HREF = "unfix.html">unfix</A>
 </P>
 <P>Computes:
 </P>
 <P><A HREF = "compute.html">compute</A>, <A HREF = "compute_modify.html">compute_modify</A>,
 <A HREF = "uncompute.html">uncompute</A>
 </P>
 <P>Output:
 </P>
 <P><A HREF = "dump.html">dump</A>, <A HREF = "dump_modify.html">dump_modify</A>,
 <A HREF = "restart.html">restart</A>, <A HREF = "thermo.html">thermo</A>,
 <A HREF = "thermo_modify.html">thermo_modify</A>, <A HREF = "thermo_style.html">thermo_style</A>,
 <A HREF = "undump.html">undump</A>, <A HREF = "write_restart.html">write_restart</A>
 </P>
 <P>Actions:
 </P>
 <P><A HREF = "delete_atoms.html">delete_atoms</A>, <A HREF = "delete_bonds.html">delete_bonds</A>,
 <A HREF = "displace_atoms.html">displace_atoms</A>,
 <A HREF = "displace_box.html">displace_box</A>, <A HREF = "minimize.html">minimize</A>,
 <A HREF = "neb.html">neb</A> <A HREF = "prd.html">prd</A>, <A HREF = "run.html">run</A>, <A HREF = "temper.html">temper</A>
 </P>
 <P>Miscellaneous:
 </P>
 <P><A HREF = "clear.html">clear</A>, <A HREF = "echo.html">echo</A>, <A HREF = "if.html">if</A>,
 <A HREF = "include.html">include</A>, <A HREF = "jump.html">jump</A>, <A HREF = "label.html">label</A>,
 <A HREF = "log.html">log</A>, <A HREF = "next.html">next</A>, <A HREF = "print.html">print</A>,
 <A HREF = "shell.html">shell</A>, <A HREF = "variable.html">variable</A>
 </P>
 <HR>
 
 <H4><A NAME = "3_5"></A><A NAME = "comm"></A>3.5 Individual commands 
 </H4>
 <P>This section lists all LAMMPS commands alphabetically, with a separate
 listing below of styles within certain commands.  The <A HREF = "#3_4">previous
 section</A> lists the same commands, grouped by category.  Note that
 some style options for some commands are part of specific LAMMPS
 packages, which means they cannot be used unless the package was
 included when LAMMPS was built.  Not all packages are included in a
 default LAMMPS build.  These dependencies are listed as Restrictions
 in the command's documentation.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "angle_coeff.html">angle_coeff</A></TD><TD ><A HREF = "angle_style.html">angle_style</A></TD><TD ><A HREF = "atom_modify.html">atom_modify</A></TD><TD ><A HREF = "atom_style.html">atom_style</A></TD><TD ><A HREF = "bond_coeff.html">bond_coeff</A></TD><TD ><A HREF = "bond_style.html">bond_style</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "boundary.html">boundary</A></TD><TD ><A HREF = "change_box.html">change_box</A></TD><TD ><A HREF = "clear.html">clear</A></TD><TD ><A HREF = "communicate.html">communicate</A></TD><TD ><A HREF = "compute.html">compute</A></TD><TD ><A HREF = "compute_modify.html">compute_modify</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "create_atoms.html">create_atoms</A></TD><TD ><A HREF = "create_box.html">create_box</A></TD><TD ><A HREF = "delete_atoms.html">delete_atoms</A></TD><TD ><A HREF = "delete_bonds.html">delete_bonds</A></TD><TD ><A HREF = "dielectric.html">dielectric</A></TD><TD ><A HREF = "dihedral_coeff.html">dihedral_coeff</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "dihedral_style.html">dihedral_style</A></TD><TD ><A HREF = "dimension.html">dimension</A></TD><TD ><A HREF = "dipole.html">dipole</A></TD><TD ><A HREF = "displace_atoms.html">displace_atoms</A></TD><TD ><A HREF = "displace_box.html">displace_box</A></TD><TD ><A HREF = "dump.html">dump</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "dump_modify.html">dump_modify</A></TD><TD ><A HREF = "echo.html">echo</A></TD><TD ><A HREF = "fix.html">fix</A></TD><TD ><A HREF = "fix_modify.html">fix_modify</A></TD><TD ><A HREF = "group.html">group</A></TD><TD ><A HREF = "if.html">if</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "improper_coeff.html">improper_coeff</A></TD><TD ><A HREF = "improper_style.html">improper_style</A></TD><TD ><A HREF = "include.html">include</A></TD><TD ><A HREF = "jump.html">jump</A></TD><TD ><A HREF = "kspace_modify.html">kspace_modify</A></TD><TD ><A HREF = "kspace_style.html">kspace_style</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "label.html">label</A></TD><TD ><A HREF = "lattice.html">lattice</A></TD><TD ><A HREF = "log.html">log</A></TD><TD ><A HREF = "mass.html">mass</A></TD><TD ><A HREF = "minimize.html">minimize</A></TD><TD ><A HREF = "min_modify.html">min_modify</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "min_style.html">min_style</A></TD><TD ><A HREF = "neb.html">neb</A></TD><TD ><A HREF = "neigh_modify.html">neigh_modify</A></TD><TD ><A HREF = "neighbor.html">neighbor</A></TD><TD ><A HREF = "newton.html">newton</A></TD><TD ><A HREF = "next.html">next</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_coeff.html">pair_coeff</A></TD><TD ><A HREF = "pair_modify.html">pair_modify</A></TD><TD ><A HREF = "pair_style.html">pair_style</A></TD><TD ><A HREF = "pair_write.html">pair_write</A></TD><TD ><A HREF = "prd.html">prd</A></TD><TD ><A HREF = "print.html">print</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "processors.html">processors</A></TD><TD ><A HREF = "read_data.html">read_data</A></TD><TD ><A HREF = "read_restart.html">read_restart</A></TD><TD ><A HREF = "region.html">region</A></TD><TD ><A HREF = "replicate.html">replicate</A></TD><TD ><A HREF = "reset_timestep.html">reset_timestep</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "restart.html">restart</A></TD><TD ><A HREF = "run.html">run</A></TD><TD ><A HREF = "run_style.html">run_style</A></TD><TD ><A HREF = "set.html">set</A></TD><TD ><A HREF = "shape.html">shape</A></TD><TD ><A HREF = "shell.html">shell</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "special_bonds.html">special_bonds</A></TD><TD ><A HREF = "temper.html">temper</A></TD><TD ><A HREF = "thermo.html">thermo</A></TD><TD ><A HREF = "thermo_modify.html">thermo_modify</A></TD><TD ><A HREF = "thermo_style.html">thermo_style</A></TD><TD ><A HREF = "timestep.html">timestep</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "uncompute.html">uncompute</A></TD><TD ><A HREF = "undump.html">undump</A></TD><TD ><A HREF = "unfix.html">unfix</A></TD><TD ><A HREF = "units.html">units</A></TD><TD ><A HREF = "variable.html">variable</A></TD><TD ><A HREF = "velocity.html">velocity</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "write_restart.html">write_restart</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Fix styles 
 </H4>
 <P>See the <A HREF = "fix.html">fix</A> command for one-line descriptions
 of each style or click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "fix_adapt.html">adapt</A></TD><TD ><A HREF = "fix_addforce.html">addforce</A></TD><TD ><A HREF = "fix_aveforce.html">aveforce</A></TD><TD ><A HREF = "fix_ave_atom.html">ave/atom</A></TD><TD ><A HREF = "fix_ave_correlate.html">ave/correlate</A></TD><TD ><A HREF = "fix_ave_histo.html">ave/histo</A></TD><TD ><A HREF = "fix_ave_spatial.html">ave/spatial</A></TD><TD ><A HREF = "fix_ave_time.html">ave/time</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_bond_break.html">bond/break</A></TD><TD ><A HREF = "fix_bond_create.html">bond/create</A></TD><TD ><A HREF = "fix_bond_swap.html">bond/swap</A></TD><TD ><A HREF = "fix_box_relax.html">box/relax</A></TD><TD ><A HREF = "fix_deform.html">deform</A></TD><TD ><A HREF = "fix_deposit.html">deposit</A></TD><TD ><A HREF = "fix_drag.html">drag</A></TD><TD ><A HREF = "fix_dt_reset.html">dt/reset</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_efield.html">efield</A></TD><TD ><A HREF = "fix_enforce2d.html">enforce2d</A></TD><TD ><A HREF = "fix_evaporate.html">evaporate</A></TD><TD ><A HREF = "fix_external.html">external</A></TD><TD ><A HREF = "fix_freeze.html">freeze</A></TD><TD ><A HREF = "fix_gravity.html">gravity</A></TD><TD ><A HREF = "fix_heat.html">heat</A></TD><TD ><A HREF = "fix_indent.html">indent</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_langevin.html">langevin</A></TD><TD ><A HREF = "fix_lineforce.html">lineforce</A></TD><TD ><A HREF = "fix_momentum.html">momentum</A></TD><TD ><A HREF = "fix_move.html">move</A></TD><TD ><A HREF = "fix_msst.html">msst</A></TD><TD ><A HREF = "fix_neb.html">neb</A></TD><TD ><A HREF = "fix_nh.html">nph</A></TD><TD ><A HREF = "fix_nph_asphere.html">nph/asphere</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_nph_sphere.html">nph/sphere</A></TD><TD ><A HREF = "fix_nh.html">npt</A></TD><TD ><A HREF = "fix_npt_asphere.html">npt/asphere</A></TD><TD ><A HREF = "fix_npt_sphere.html">npt/sphere</A></TD><TD ><A HREF = "fix_nve.html">nve</A></TD><TD ><A HREF = "fix_nve_asphere.html">nve/asphere</A></TD><TD ><A HREF = "fix_nve_limit.html">nve/limit</A></TD><TD ><A HREF = "fix_nve_noforce.html">nve/noforce</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_nve_sphere.html">nve/sphere</A></TD><TD ><A HREF = "fix_nh.html">nvt</A></TD><TD ><A HREF = "fix_nvt_asphere.html">nvt/asphere</A></TD><TD ><A HREF = "fix_nvt_sllod.html">nvt/sllod</A></TD><TD ><A HREF = "fix_nvt_sphere.html">nvt/sphere</A></TD><TD ><A HREF = "fix_orient_fcc.html">orient/fcc</A></TD><TD ><A HREF = "fix_planeforce.html">planeforce</A></TD><TD ><A HREF = "fix_poems.html">poems</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_pour.html">pour</A></TD><TD ><A HREF = "fix_press_berendsen.html">press/berendsen</A></TD><TD ><A HREF = "fix_print.html">print</A></TD><TD ><A HREF = "fix_qeq_comb.html">qeq/comb</A></TD><TD ><A HREF = "fix_reax_bonds.html">reax/bonds</A></TD><TD ><A HREF = "fix_recenter.html">recenter</A></TD><TD ><A HREF = "fix_rigid.html">rigid</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nve</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_rigid.html">rigid/nvt</A></TD><TD ><A HREF = "fix_setforce.html">setforce</A></TD><TD ><A HREF = "fix_shake.html">shake</A></TD><TD ><A HREF = "fix_spring.html">spring</A></TD><TD ><A HREF = "fix_spring_rg.html">spring/rg</A></TD><TD ><A HREF = "fix_spring_self.html">spring/self</A></TD><TD ><A HREF = "fix_srd.html">srd</A></TD><TD ><A HREF = "fix_store_force.html">store/force</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_store_state.html">store/state</A></TD><TD ><A HREF = "fix_temp_berendsen.html">temp/berendsen</A></TD><TD ><A HREF = "fix_temp_rescale.html">temp/rescale</A></TD><TD ><A HREF = "fix_thermal_conductivity.html">thermal/conductivity</A></TD><TD ><A HREF = "fix_tmd.html">tmd</A></TD><TD ><A HREF = "fix_ttm.html">ttm</A></TD><TD ><A HREF = "fix_viscosity.html">viscosity</A></TD><TD ><A HREF = "fix_viscous.html">viscous</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_wall.html">wall/colloid</A></TD><TD ><A HREF = "fix_wall_gran.html">wall/gran</A></TD><TD ><A HREF = "fix_wall.html">wall/harmonic</A></TD><TD ><A HREF = "fix_wall.html">wall/lj126</A></TD><TD ><A HREF = "fix_wall.html">wall/lj93</A></TD><TD ><A HREF = "fix_wall_reflect.html">wall/reflect</A></TD><TD ><A HREF = "fix_wall_region.html">wall/region</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are fix styles contributed by users, which can be used if
 <A HREF = "Section_start.html#2_3">LAMMPS is built with the appropriate package</A>.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "fix_atc.html">atc</A></TD><TD ><A HREF = "fix_imd.html">imd</A></TD><TD ><A HREF = "fix_langevin_eff.html">langevin/eff</A></TD><TD ><A HREF = "fix_nh_eff.html">nph/eff</A></TD><TD ><A HREF = "fix_nh_eff.html">npt/eff</A></TD><TD ><A HREF = "fix_nve_eff.html">nve/eff</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "fix_nh_eff.html">nvt/eff</A></TD><TD ><A HREF = "fix_nvt_sllod_eff.html">nvt/sllod/eff</A></TD><TD ><A HREF = "fix_qeq_reax.html">qeq/reax</A></TD><TD ><A HREF = "fix_smd.html">smd</A></TD><TD ><A HREF = "fix_temp_rescale_eff.html">temp/rescale/eff</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Compute styles 
 </H4>
 <P>See the <A HREF = "compute.html">compute</A> command for one-line descriptions of
 each style or click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "compute_angle_local.html">angle/local</A></TD><TD ><A HREF = "compute_atom_molecule.html">atom/molecule</A></TD><TD ><A HREF = "compute_bond_local.html">bond/local</A></TD><TD ><A HREF = "compute_centro_atom.html">centro/atom</A></TD><TD ><A HREF = "compute_cna_atom.html">cna/atom</A></TD><TD ><A HREF = "compute_com.html">com</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_com_molecule.html">com/molecule</A></TD><TD ><A HREF = "compute_coord_atom.html">coord/atom</A></TD><TD ><A HREF = "compute_damage_atom.html">damage/atom</A></TD><TD ><A HREF = "compute_dihedral_local.html">dihedral/local</A></TD><TD ><A HREF = "compute_displace_atom.html">displace/atom</A></TD><TD ><A HREF = "compute_erotate_asphere.html">erotate/asphere</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_erotate_sphere.html">erotate/sphere</A></TD><TD ><A HREF = "compute_event_displace.html">event/displace</A></TD><TD ><A HREF = "compute_group_group.html">group/group</A></TD><TD ><A HREF = "compute_gyration.html">gyration</A></TD><TD ><A HREF = "compute_gyration_molecule.html">gyration/molecule</A></TD><TD ><A HREF = "compute_heat_flux.html">heat/flux</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_improper_local.html">improper/local</A></TD><TD ><A HREF = "compute_ke.html">ke</A></TD><TD ><A HREF = "compute_ke_atom.html">ke/atom</A></TD><TD ><A HREF = "compute_msd.html">msd</A></TD><TD ><A HREF = "compute_msd_molecule.html">msd/molecule</A></TD><TD ><A HREF = "compute_pair.html">pair</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_pair_local.html">pair/local</A></TD><TD ><A HREF = "compute_pe.html">pe</A></TD><TD ><A HREF = "compute_pe_atom.html">pe/atom</A></TD><TD ><A HREF = "compute_pressure.html">pressure</A></TD><TD ><A HREF = "compute_property_atom.html">property/atom</A></TD><TD ><A HREF = "compute_property_local.html">property/local</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_property_molecule.html">property/molecule</A></TD><TD ><A HREF = "compute_rdf.html">rdf</A></TD><TD ><A HREF = "compute_reduce.html">reduce</A></TD><TD ><A HREF = "compute_reduce.html">reduce/region</A></TD><TD ><A HREF = "compute_stress_atom.html">stress/atom</A></TD><TD ><A HREF = "compute_temp.html">temp</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_temp_asphere.html">temp/asphere</A></TD><TD ><A HREF = "compute_temp_com.html">temp/com</A></TD><TD ><A HREF = "compute_temp_deform.html">temp/deform</A></TD><TD ><A HREF = "compute_temp_partial.html">temp/partial</A></TD><TD ><A HREF = "compute_temp_profile.html">temp/profile</A></TD><TD ><A HREF = "compute_temp_ramp.html">temp/ramp</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_temp_region.html">temp/region</A></TD><TD ><A HREF = "compute_temp_sphere.html">temp/sphere</A></TD><TD ><A HREF = "compute_ti.html">ti</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are compute styles contributed by users, which can be used if
 <A HREF = "Section_start.html#2_3">LAMMPS is built with the appropriate package</A>.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "compute_ackland_atom.html">ackland/atom</A></TD><TD ><A HREF = "compute_ke_eff.html">ke/eff</A></TD><TD ><A HREF = "compute_ke_atom_eff.html">ke/atom/eff</A></TD><TD ><A HREF = "compute_temp_eff.html">temp/eff</A></TD><TD ><A HREF = "compute_temp_deform_eff.html">temp/deform/eff</A></TD><TD ><A HREF = "compute_temp_region_eff.html">temp/region/eff</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Pair_style potentials 
 </H4>
 <P>See the <A HREF = "pair_style.html">pair_style</A> command for an overview of pair
 potentials.  Click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "pair_none.html">none</A></TD><TD ><A HREF = "pair_hybrid.html">hybrid</A></TD><TD ><A HREF = "pair_hybrid.html">hybrid/overlay</A></TD><TD ><A HREF = "pair_airebo.html">airebo</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_born.html">born</A></TD><TD ><A HREF = "pair_born.html">born/coul/long</A></TD><TD ><A HREF = "pair_buck.html">buck</A></TD><TD ><A HREF = "pair_buck.html">buck/coul/cut</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_buck.html">buck/coul/long</A></TD><TD ><A HREF = "pair_colloid.html">colloid</A></TD><TD ><A HREF = "pair_comb.html">comb</A></TD><TD ><A HREF = "pair_coul.html">coul/cut</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_coul.html">coul/debye</A></TD><TD ><A HREF = "pair_coul.html">coul/long</A></TD><TD ><A HREF = "pair_dipole.html">dipole/cut</A></TD><TD ><A HREF = "pair_dpd.html">dpd</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_dpd.html">dpd/tstat</A></TD><TD ><A HREF = "pair_dsmc.html">dsmc</A></TD><TD ><A HREF = "pair_eam.html">eam</A></TD><TD ><A HREF = "pair_eam.html">eam/opt</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_eam.html">eam/alloy</A></TD><TD ><A HREF = "pair_eam.html">eam/alloy/opt</A></TD><TD ><A HREF = "pair_eam.html">eam/fs</A></TD><TD ><A HREF = "pair_eam.html">eam/fs/opt</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_eim.html">eim</A></TD><TD ><A HREF = "pair_gauss.html">gauss</A></TD><TD ><A HREF = "pair_gayberne.html">gayberne</A></TD><TD ><A HREF = "pair_gayberne.html">gayberne/gpu</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_gran.html">gran/hertz/history</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke/history</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm/implicit</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long/opt</A></TD><TD ><A HREF = "pair_class2.html">lj/class2</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_class2.html">lj/class2/coul/cut</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/gpu</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/opt</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/debye</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/tip4p</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_lj_smooth.html">lj/smooth</A></TD><TD ><A HREF = "pair_lj96_cut.html">lj96/cut</A></TD><TD ><A HREF = "pair_lubricate.html">lubricate</A></TD><TD ><A HREF = "pair_meam.html">meam</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_morse.html">morse</A></TD><TD ><A HREF = "pair_morse.html">morse/opt</A></TD><TD ><A HREF = "pair_peri.html">peri/lps</A></TD><TD ><A HREF = "pair_peri.html">peri/pmb</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_resquared.html">resquared</A></TD><TD ><A HREF = "pair_soft.html">soft</A></TD><TD ><A HREF = "pair_sw.html">sw</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_table.html">table</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff</A></TD><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid</A> 
+<TR ALIGN="center"><TD ><A HREF = "pair_gran.html">gran/hertz/history</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke/history</A></TD><TD ><A HREF = "pair_hbond_dreiding.html">hbond/dreiding/lj</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_hbond_dreiding.html">hbond/dreiding/morse</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm/implicit</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long/opt</A></TD><TD ><A HREF = "pair_class2.html">lj/class2</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/cut</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/long</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/gpu</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/opt</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/debye</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/tip4p</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_gromacs.html">lj/gromacs</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs</A></TD><TD ><A HREF = "pair_lj_smooth.html">lj/smooth</A></TD><TD ><A HREF = "pair_lj96_cut.html">lj96/cut</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_lubricate.html">lubricate</A></TD><TD ><A HREF = "pair_meam.html">meam</A></TD><TD ><A HREF = "pair_morse.html">morse</A></TD><TD ><A HREF = "pair_morse.html">morse/opt</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_peri.html">peri/lps</A></TD><TD ><A HREF = "pair_peri.html">peri/pmb</A></TD><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_resquared.html">resquared</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_soft.html">soft</A></TD><TD ><A HREF = "pair_sw.html">sw</A></TD><TD ><A HREF = "pair_table.html">table</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa</A></TD><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are pair styles contributed by users, which can be used if
 <A HREF = "Section_start.html#2_3">LAMMPS is built with the appropriate package</A>.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "pair_buck_coul.html">buck/coul</A></TD><TD ><A HREF = "pair_cmm.html">cg/cmm</A></TD><TD ><A HREF = "pair_cmm.html">cg/cmm/coul/cut</A></TD><TD ><A HREF = "pair_cmm.html">cg/cmm/coul/long</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_eam.html">eam/cd</A></TD><TD ><A HREF = "pair_eff.html">eff/cut</A></TD><TD ><A HREF = "pair_lj_coul.html">lj/coul</A></TD><TD ><A HREF = "pair_reax_c.html">reax/c</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Bond_style potentials 
 </H4>
 <P>See the <A HREF = "bond_style.html">bond_style</A> command for an overview of bond
 potentials.  Click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "bond_none.html">none</A></TD><TD WIDTH="100"><A HREF = "bond_hybrid.html">hybrid</A></TD><TD WIDTH="100"><A HREF = "bond_class2.html">class2</A></TD><TD WIDTH="100"><A HREF = "bond_fene.html">fene</A></TD></TR>
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "bond_fene_expand.html">fene/expand</A></TD><TD WIDTH="100"><A HREF = "bond_harmonic.html">harmonic</A></TD><TD WIDTH="100"><A HREF = "bond_morse.html">morse</A></TD><TD WIDTH="100"><A HREF = "bond_nonlinear.html">nonlinear</A></TD></TR>
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "bond_quartic.html">quartic</A></TD><TD WIDTH="100"><A HREF = "bond_table.html">table</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Angle_style potentials 
 </H4>
 <P>See the <A HREF = "angle_style.html">angle_style</A> command for an overview of
 angle potentials.  Click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "angle_none.html">none</A></TD><TD WIDTH="100"><A HREF = "angle_hybrid.html">hybrid</A></TD><TD WIDTH="100"><A HREF = "angle_charmm.html">charmm</A></TD><TD WIDTH="100"><A HREF = "angle_class2.html">class2</A></TD></TR>
-<TR ALIGN="center"><TD WIDTH="100"><A HREF = "angle_cosine.html">cosine</A></TD><TD WIDTH="100"><A HREF = "angle_cosine_delta.html">cosine/delta</A></TD><TD WIDTH="100"><A HREF = "angle_cosine_squared.html">cosine/squared</A></TD><TD WIDTH="100"><A HREF = "angle_harmonic.html">harmonic</A></TD></TR>
-<TR ALIGN="center"><TD WIDTH="100"><A HREF = "angle_table.html">table</A> 
+<TR ALIGN="center"><TD WIDTH="100"><A HREF = "angle_cosine.html">cosine</A></TD><TD WIDTH="100"><A HREF = "angle_cosine_delta.html">cosine/delta</A></TD><TD WIDTH="100"><A HREF = "angle_cosine_periodic.html">cosine/periodic</A></TD><TD WIDTH="100"><A HREF = "angle_cosine_squared.html">cosine/squared</A></TD></TR>
+<TR ALIGN="center"><TD WIDTH="100"><A HREF = "angle_harmonic.html">harmonic</A></TD><TD WIDTH="100"><A HREF = "angle_table.html">table</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are angle styles contributed by users, which can be used if
 <A HREF = "Section_start.html#2_3">LAMMPS is built with the appropriate package</A>.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD ><A HREF = "angle_cmm.html">cg/cmm</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Dihedral_style potentials 
 </H4>
 <P>See the <A HREF = "dihedral_style.html">dihedral_style</A> command for an overview
 of dihedral potentials.  Click on the style itself for a full
 description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "dihedral_none.html">none</A></TD><TD WIDTH="100"><A HREF = "dihedral_hybrid.html">hybrid</A></TD><TD WIDTH="100"><A HREF = "dihedral_charmm.html">charmm</A></TD><TD WIDTH="100"><A HREF = "dihedral_class2.html">class2</A></TD></TR>
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "dihedral_harmonic.html">harmonic</A></TD><TD WIDTH="100"><A HREF = "dihedral_helix.html">helix</A></TD><TD WIDTH="100"><A HREF = "dihedral_multi_harmonic.html">multi/harmonic</A></TD><TD WIDTH="100"><A HREF = "dihedral_opls.html">opls</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Improper_style potentials 
 </H4>
 <P>See the <A HREF = "improper_style.html">improper_style</A> command for an overview
 of improper potentials.  Click on the style itself for a full
 description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "improper_none.html">none</A></TD><TD WIDTH="100"><A HREF = "improper_hybrid.html">hybrid</A></TD><TD WIDTH="100"><A HREF = "improper_class2.html">class2</A></TD><TD WIDTH="100"><A HREF = "improper_cvff.html">cvff</A></TD></TR>
-<TR ALIGN="center"><TD WIDTH="100"><A HREF = "improper_harmonic.html">harmonic</A> 
+<TR ALIGN="center"><TD WIDTH="100"><A HREF = "improper_harmonic.html">harmonic</A></TD><TD WIDTH="100"><A HREF = "improper_umbrella.html">umbrella</A> 
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <H4>Kspace solvers 
 </H4>
 <P>See the <A HREF = "kspace_style.html">kspace_style</A> command for an overview of
 Kspace solvers.  Click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "kspace_style.html">ewald</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm/tip4p</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are Kspace solvers contributed by users, which can be used if
 <A HREF = "Section_start.html#2_3">LAMMPS is built with the appropriate package</A>.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR ALIGN="center"><TD WIDTH="100"><A HREF = "kspace_style.html">ewald/n</A> 
 </TD></TR></TABLE></DIV>
 
 </HTML>
diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt
index 9dfb33cff..53997e819 100644
--- a/doc/Section_commands.txt
+++ b/doc/Section_commands.txt
@@ -1,714 +1,718 @@
 "Previous Section"_Section_start.html - "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc - "Next Section"_Section_howto.html :c
 
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 
 :line
 
 3. Commands :h3
 
 This section describes how a LAMMPS input script is formatted and what
 commands are used to define a LAMMPS simulation.
 
 3.1 "LAMMPS input script"_#3_1
 3.2 "Parsing rules"_#3_2
 3.3 "Input script structure"_#3_3
 3.4 "Commands listed by category"_#3_4
 3.5 "Commands listed alphabetically"_#3_5 :all(b)
 
 :line
 
 3.1 LAMMPS input script :link(3_1),h4
 
 LAMMPS executes by reading commands from a input script (text file),
 one line at a time.  When the input script ends, LAMMPS exits.  Each
 command causes LAMMPS to take some action.  It may set an internal
 variable, read in a file, or run a simulation.  Most commands have
 default settings, which means you only need to use the command if you
 wish to change the default.
 
 In many cases, the ordering of commands in an input script is not
 important.  However the following rules apply:
 
 (1) LAMMPS does not read your entire input script and then perform a
 simulation with all the settings.  Rather, the input script is read
 one line at a time and each command takes effect when it is read.
 Thus this sequence of commands:
 
 timestep 0.5 
 run      100 
 run      100 :pre
 
 does something different than this sequence:
 
 run      100 
 timestep 0.5 
 run      100 :pre
 
 In the first case, the specified timestep (0.5 fmsec) is used for two
 simulations of 100 timesteps each.  In the 2nd case, the default
 timestep (1.0 fmsec) is used for the 1st 100 step simulation and a 0.5
 fmsec timestep is used for the 2nd one.
 
 (2) Some commands are only valid when they follow other commands.  For
 example you cannot set the temperature of a group of atoms until atoms
 have been defined and a group command is used to define which atoms
 belong to the group.
 
 (3) Sometimes command B will use values that can be set by command A.
 This means command A must precede command B in the input script if it
 is to have the desired effect.  For example, the
 "read_data"_read_data.html command initializes the system by setting
 up the simulation box and assigning atoms to processors.  If default
 values are not desired, the "processors"_processors.html and
 "boundary"_boundary.html commands need to be used before read_data to
 tell LAMMPS how to map processors to the simulation box.
 
 Many input script errors are detected by LAMMPS and an ERROR or
 WARNING message is printed.  "This section"_Section_errors.html gives
 more information on what errors mean.  The documentation for each
 command lists restrictions on how the command can be used.
 
 :line
 
 3.2 Parsing rules :link(3_2),h4
 
 Each non-blank line in the input script is treated as a command.
 LAMMPS commands are case sensitive.  Command names are lower-case, as
 are specified command arguments.  Upper case letters may be used in
 file names or user-chosen ID strings.
 
 Here is how each line in the input script is parsed by LAMMPS:
 
 (1) If the last printable character on the line is a "&" character
 (with no surrounding quotes), the command is assumed to continue on
 the next line.  The next line is concatenated to the previous line by
 removing the "&" character and newline.  This allows long commands to
 be continued across two or more lines.
 
 (2) All characters from the first "#" character onward are treated as
 comment and discarded.  See an exception in (6).  Note that a
 comment after a trailing "&" character will prevent the command from
 continuing on the next line.  Also note that for multi-line commands a
 single leading "#" will comment out the entire command.
 
 (3) The line is searched repeatedly for $ characters, which indicate
 variables that are replaced with a text string.  See an exception in
 (6).  If the $ is followed by curly brackets, then the variable name
 is the text inside the curly brackets.  If no curly brackets follow
 the $, then the variable name is the single character immediately
 following the $.  Thus $\{myTemp\} and $x refer to variable names
 "myTemp" and "x".  See the "variable"_variable.html command for
 details of how strings are assigned to variables and how they are
 substituted for in input script commands.
 
 (4) The line is broken into "words" separated by whitespace (tabs,
 spaces).  Note that words can thus contain letters, digits,
 underscores, or punctuation characters.
 
 (5) The first word is the command name.  All successive words in the
 line are arguments.
 
 (6) If you want text with spaces to be treated as a single argument,
 it can be enclosed in double quotes.  The quotes are removed when the
 single argument is stored in its final form.  See the "dump modify
 format"_dump_modify.html or "if"_if.html commands for examples.  A "#"
 or "$" character that is between double quotes will not be treated as
 a comment indicator in (2) or substituted for as a variable in (3).
 
 :line
 
 3.3 Input script structure :h4,link(3_3)
 
 This section describes the structure of a typical LAMMPS input script.
 The "examples" directory in the LAMMPS distribution contains many
 sample input scripts; the corresponding problems are discussed in
 "this section"_Section_example.html, and animated on the "LAMMPS WWW
 Site"_lws.
 
 A LAMMPS input script typically has 4 parts:
 
 Initialization
 Atom definition
 Settings
 Run a simulation :ol
 
 The last 2 parts can be repeated as many times as desired.  I.e. run a
 simulation, change some settings, run some more, etc.  Each of the 4
 parts is now described in more detail.  Remember that almost all the
 commands need only be used if a non-default value is desired.
 
 (1) Initialization
 
 Set parameters that need to be defined before atoms are created or
 read-in from a file.
 
 The relevant commands are "units"_units.html,
 "dimension"_dimension.html, "newton"_newton.html,
 "processors"_processors.html, "boundary"_boundary.html,
 "atom_style"_atom_style.html, "atom_modify"_atom_modify.html.
 
 If force-field parameters appear in the files that will be read, these
 commands tell LAMMPS what kinds of force fields are being used:
 "pair_style"_pair_style.html, "bond_style"_bond_style.html,
 "angle_style"_angle_style.html, "dihedral_style"_dihedral_style.html,
 "improper_style"_improper_style.html.
 
 (2) Atom definition
 
 There are 3 ways to define atoms in LAMMPS.  Read them in from a data
 or restart file via the "read_data"_read_data.html or
 "read_restart"_read_restart.html commands.  These files can contain
 molecular topology information.  Or create atoms on a lattice (with no
 molecular topology), using these commands: "lattice"_lattice.html,
 "region"_region.html, "create_box"_create_box.html,
 "create_atoms"_create_atoms.html.  The entire set of atoms can be
 duplicated to make a larger simulation using the
 "replicate"_replicate.html command.
 
 (3) Settings
 
 Once atoms and molecular topology are defined, a variety of settings
 can be specified: force field coefficients, simulation parameters,
 output options, etc.
 
 Force field coefficients are set by these commands (they can also be
 set in the read-in files): "pair_coeff"_pair_coeff.html,
 "bond_coeff"_bond_coeff.html, "angle_coeff"_angle_coeff.html,
 "dihedral_coeff"_dihedral_coeff.html,
 "improper_coeff"_improper_coeff.html,
 "kspace_style"_kspace_style.html, "dielectric"_dielectric.html,
 "special_bonds"_special_bonds.html.
 
 Various simulation parameters are set by these commands:
 "neighbor"_neighbor.html, "neigh_modify"_neigh_modify.html,
 "group"_group.html, "timestep"_timestep.html,
 "reset_timestep"_reset_timestep.html, "run_style"_run_style.html,
 "min_style"_min_style.html, "min_modify"_min_modify.html.
 
 Fixes impose a variety of boundary conditions, time integration, and
 diagnostic options.  The "fix"_fix.html command comes in many flavors.
 
 Various computations can be specified for execution during a
 simulation using the "compute"_compute.html,
 "compute_modify"_compute_modify.html, and "variable"_variable.html
 commands.
 
 Output options are set by the "thermo"_thermo.html, "dump"_dump.html,
 and "restart"_restart.html commands.
 
 (4) Run a simulation
 
 A molecular dynamics simulation is run using the "run"_run.html
 command.  Energy minimization (molecular statics) is performed using
 the "minimize"_minimize.html command.  A parallel tempering
 (replica-exchange) simulation can be run using the
 "temper"_temper.html command.
 
 :line
 
 3.4 Commands listed by category :link(3_4),h4
 
 This section lists all LAMMPS commands, grouped by category.  The
 "next section"_#3_5 lists the same commands alphabetically.  Note that
 some style options for some commands are part of specific LAMMPS
 packages, which means they cannot be used unless the package was
 included when LAMMPS was built.  Not all packages are included in a
 default LAMMPS build.  These dependencies are listed as Restrictions
 in the command's documentation.
 
 Initialization:
 
 "atom_modify"_atom_modify.html, "atom_style"_atom_style.html,
 "boundary"_boundary.html, "dimension"_dimension.html,
 "newton"_newton.html, "processors"_processors.html, "units"_units.html
 
 Atom definition:
 
 "create_atoms"_create_atoms.html, "create_box"_create_box.html,
 "lattice"_lattice.html, "read_data"_read_data.html,
 "read_restart"_read_restart.html, "region"_region.html,
 "replicate"_replicate.html
 
 Force fields:
 
 "angle_coeff"_angle_coeff.html, "angle_style"_angle_style.html,
 "bond_coeff"_bond_coeff.html, "bond_style"_bond_style.html,
 "dielectric"_dielectric.html, "dihedral_coeff"_dihedral_coeff.html,
 "dihedral_style"_dihedral_style.html,
 "improper_coeff"_improper_coeff.html,
 "improper_style"_improper_style.html,
 "kspace_modify"_kspace_modify.html, "kspace_style"_kspace_style.html,
 "pair_coeff"_pair_coeff.html, "pair_modify"_pair_modify.html,
 "pair_style"_pair_style.html, "pair_write"_pair_write.html,
 "special_bonds"_special_bonds.html
 
 Settings:
 
 "communicate"_communicate.html, "dipole"_dipole.html,
 "group"_group.html, "mass"_mass.html, "min_modify"_min_modify.html,
 "min_style"_min_style.html, "neigh_modify"_neigh_modify.html,
 "neighbor"_neighbor.html, "reset_timestep"_reset_timestep.html,
 "run_style"_run_style.html, "set"_set.html, "shape"_shape.html,
 "timestep"_timestep.html, "velocity"_velocity.html
 
 Fixes:
 
 "fix"_fix.html, "fix_modify"_fix_modify.html, "unfix"_unfix.html
 
 Computes:
 
 "compute"_compute.html, "compute_modify"_compute_modify.html,
 "uncompute"_uncompute.html
 
 Output:
 
 "dump"_dump.html, "dump_modify"_dump_modify.html,
 "restart"_restart.html, "thermo"_thermo.html,
 "thermo_modify"_thermo_modify.html, "thermo_style"_thermo_style.html,
 "undump"_undump.html, "write_restart"_write_restart.html
 
 Actions:
 
 "delete_atoms"_delete_atoms.html, "delete_bonds"_delete_bonds.html,
 "displace_atoms"_displace_atoms.html,
 "displace_box"_displace_box.html, "minimize"_minimize.html,
 "neb"_neb.html "prd"_prd.html, "run"_run.html, "temper"_temper.html
 
 Miscellaneous:
 
 "clear"_clear.html, "echo"_echo.html, "if"_if.html,
 "include"_include.html, "jump"_jump.html, "label"_label.html,
 "log"_log.html, "next"_next.html, "print"_print.html,
 "shell"_shell.html, "variable"_variable.html
 
 :line
 
 3.5 Individual commands :h4,link(3_5),link(comm)
 
 This section lists all LAMMPS commands alphabetically, with a separate
 listing below of styles within certain commands.  The "previous
 section"_#3_4 lists the same commands, grouped by category.  Note that
 some style options for some commands are part of specific LAMMPS
 packages, which means they cannot be used unless the package was
 included when LAMMPS was built.  Not all packages are included in a
 default LAMMPS build.  These dependencies are listed as Restrictions
 in the command's documentation.
 
 "angle_coeff"_angle_coeff.html,
 "angle_style"_angle_style.html,
 "atom_modify"_atom_modify.html,
 "atom_style"_atom_style.html,
 "bond_coeff"_bond_coeff.html,
 "bond_style"_bond_style.html,
 "boundary"_boundary.html,
 "change_box"_change_box.html,
 "clear"_clear.html,
 "communicate"_communicate.html,
 "compute"_compute.html,
 "compute_modify"_compute_modify.html,
 "create_atoms"_create_atoms.html,
 "create_box"_create_box.html,
 "delete_atoms"_delete_atoms.html,
 "delete_bonds"_delete_bonds.html,
 "dielectric"_dielectric.html,
 "dihedral_coeff"_dihedral_coeff.html,
 "dihedral_style"_dihedral_style.html,
 "dimension"_dimension.html,
 "dipole"_dipole.html,
 "displace_atoms"_displace_atoms.html,
 "displace_box"_displace_box.html,
 "dump"_dump.html,
 "dump_modify"_dump_modify.html,
 "echo"_echo.html,
 "fix"_fix.html,
 "fix_modify"_fix_modify.html,
 "group"_group.html,
 "if"_if.html,
 "improper_coeff"_improper_coeff.html,
 "improper_style"_improper_style.html,
 "include"_include.html,
 "jump"_jump.html,
 "kspace_modify"_kspace_modify.html,
 "kspace_style"_kspace_style.html,
 "label"_label.html,
 "lattice"_lattice.html,
 "log"_log.html,
 "mass"_mass.html,
 "minimize"_minimize.html,
 "min_modify"_min_modify.html,
 "min_style"_min_style.html,
 "neb"_neb.html,
 "neigh_modify"_neigh_modify.html,
 "neighbor"_neighbor.html,
 "newton"_newton.html,
 "next"_next.html,
 "pair_coeff"_pair_coeff.html,
 "pair_modify"_pair_modify.html,
 "pair_style"_pair_style.html,
 "pair_write"_pair_write.html,
 "prd"_prd.html,
 "print"_print.html,
 "processors"_processors.html,
 "read_data"_read_data.html,
 "read_restart"_read_restart.html,
 "region"_region.html,
 "replicate"_replicate.html,
 "reset_timestep"_reset_timestep.html,
 "restart"_restart.html,
 "run"_run.html,
 "run_style"_run_style.html,
 "set"_set.html,
 "shape"_shape.html,
 "shell"_shell.html,
 "special_bonds"_special_bonds.html,
 "temper"_temper.html,
 "thermo"_thermo.html,
 "thermo_modify"_thermo_modify.html,
 "thermo_style"_thermo_style.html,
 "timestep"_timestep.html,
 "uncompute"_uncompute.html,
 "undump"_undump.html,
 "unfix"_unfix.html,
 "units"_units.html,
 "variable"_variable.html,
 "velocity"_velocity.html,
 "write_restart"_write_restart.html :tb(c=6,ea=c)
 
 :line
 
 Fix styles :h4
 
 See the "fix"_fix.html command for one-line descriptions
 of each style or click on the style itself for a full description:
 
 "adapt"_fix_adapt.html,
 "addforce"_fix_addforce.html,
 "aveforce"_fix_aveforce.html,
 "ave/atom"_fix_ave_atom.html,
 "ave/correlate"_fix_ave_correlate.html,
 "ave/histo"_fix_ave_histo.html,
 "ave/spatial"_fix_ave_spatial.html,
 "ave/time"_fix_ave_time.html,
 "bond/break"_fix_bond_break.html,
 "bond/create"_fix_bond_create.html,
 "bond/swap"_fix_bond_swap.html,
 "box/relax"_fix_box_relax.html,
 "deform"_fix_deform.html,
 "deposit"_fix_deposit.html,
 "drag"_fix_drag.html,
 "dt/reset"_fix_dt_reset.html,
 "efield"_fix_efield.html,
 "enforce2d"_fix_enforce2d.html,
 "evaporate"_fix_evaporate.html,
 "external"_fix_external.html,
 "freeze"_fix_freeze.html,
 "gravity"_fix_gravity.html,
 "heat"_fix_heat.html,
 "indent"_fix_indent.html,
 "langevin"_fix_langevin.html,
 "lineforce"_fix_lineforce.html,
 "momentum"_fix_momentum.html,
 "move"_fix_move.html,
 "msst"_fix_msst.html,
 "neb"_fix_neb.html,
 "nph"_fix_nh.html,
 "nph/asphere"_fix_nph_asphere.html,
 "nph/sphere"_fix_nph_sphere.html,
 "npt"_fix_nh.html,
 "npt/asphere"_fix_npt_asphere.html,
 "npt/sphere"_fix_npt_sphere.html,
 "nve"_fix_nve.html,
 "nve/asphere"_fix_nve_asphere.html,
 "nve/limit"_fix_nve_limit.html,
 "nve/noforce"_fix_nve_noforce.html,
 "nve/sphere"_fix_nve_sphere.html,
 "nvt"_fix_nh.html,
 "nvt/asphere"_fix_nvt_asphere.html,
 "nvt/sllod"_fix_nvt_sllod.html,
 "nvt/sphere"_fix_nvt_sphere.html,
 "orient/fcc"_fix_orient_fcc.html,
 "planeforce"_fix_planeforce.html,
 "poems"_fix_poems.html,
 "pour"_fix_pour.html,
 "press/berendsen"_fix_press_berendsen.html,
 "print"_fix_print.html,
 "qeq/comb"_fix_qeq_comb.html,
 "reax/bonds"_fix_reax_bonds.html,
 "recenter"_fix_recenter.html,
 "rigid"_fix_rigid.html,
 "rigid/nve"_fix_rigid.html,
 "rigid/nvt"_fix_rigid.html,
 "setforce"_fix_setforce.html,
 "shake"_fix_shake.html,
 "spring"_fix_spring.html,
 "spring/rg"_fix_spring_rg.html,
 "spring/self"_fix_spring_self.html,
 "srd"_fix_srd.html,
 "store/force"_fix_store_force.html,
 "store/state"_fix_store_state.html,
 "temp/berendsen"_fix_temp_berendsen.html,
 "temp/rescale"_fix_temp_rescale.html,
 "thermal/conductivity"_fix_thermal_conductivity.html,
 "tmd"_fix_tmd.html,
 "ttm"_fix_ttm.html,
 "viscosity"_fix_viscosity.html,
 "viscous"_fix_viscous.html,
 "wall/colloid"_fix_wall.html,
 "wall/gran"_fix_wall_gran.html,
 "wall/harmonic"_fix_wall.html,
 "wall/lj126"_fix_wall.html,
 "wall/lj93"_fix_wall.html,
 "wall/reflect"_fix_wall_reflect.html,
 "wall/region"_fix_wall_region.html :tb(c=8,ea=c)
 
 These are fix styles contributed by users, which can be used if
 "LAMMPS is built with the appropriate package"_Section_start.html#2_3.
 
 "atc"_fix_atc.html,
 "imd"_fix_imd.html,
 "langevin/eff"_fix_langevin_eff.html,
 "nph/eff"_fix_nh_eff.html,
 "npt/eff"_fix_nh_eff.html,
 "nve/eff"_fix_nve_eff.html,
 "nvt/eff"_fix_nh_eff.html,
 "nvt/sllod/eff"_fix_nvt_sllod_eff.html,
 "qeq/reax"_fix_qeq_reax.html,
 "smd"_fix_smd.html,
 "temp/rescale/eff"_fix_temp_rescale_eff.html :tb(c=6,ea=c)
 
 :line
 
 Compute styles :h4
 
 See the "compute"_compute.html command for one-line descriptions of
 each style or click on the style itself for a full description:
 
 "angle/local"_compute_angle_local.html,
 "atom/molecule"_compute_atom_molecule.html,
 "bond/local"_compute_bond_local.html,
 "centro/atom"_compute_centro_atom.html,
 "cna/atom"_compute_cna_atom.html,
 "com"_compute_com.html,
 "com/molecule"_compute_com_molecule.html,
 "coord/atom"_compute_coord_atom.html,
 "damage/atom"_compute_damage_atom.html,
 "dihedral/local"_compute_dihedral_local.html,
 "displace/atom"_compute_displace_atom.html,
 "erotate/asphere"_compute_erotate_asphere.html,
 "erotate/sphere"_compute_erotate_sphere.html,
 "event/displace"_compute_event_displace.html,
 "group/group"_compute_group_group.html,
 "gyration"_compute_gyration.html,
 "gyration/molecule"_compute_gyration_molecule.html,
 "heat/flux"_compute_heat_flux.html,
 "improper/local"_compute_improper_local.html,
 "ke"_compute_ke.html,
 "ke/atom"_compute_ke_atom.html,
 "msd"_compute_msd.html,
 "msd/molecule"_compute_msd_molecule.html,
 "pair"_compute_pair.html,
 "pair/local"_compute_pair_local.html,
 "pe"_compute_pe.html,
 "pe/atom"_compute_pe_atom.html,
 "pressure"_compute_pressure.html,
 "property/atom"_compute_property_atom.html,
 "property/local"_compute_property_local.html,
 "property/molecule"_compute_property_molecule.html,
 "rdf"_compute_rdf.html,
 "reduce"_compute_reduce.html,
 "reduce/region"_compute_reduce.html,
 "stress/atom"_compute_stress_atom.html,
 "temp"_compute_temp.html,
 "temp/asphere"_compute_temp_asphere.html,
 "temp/com"_compute_temp_com.html,
 "temp/deform"_compute_temp_deform.html,
 "temp/partial"_compute_temp_partial.html,
 "temp/profile"_compute_temp_profile.html,
 "temp/ramp"_compute_temp_ramp.html,
 "temp/region"_compute_temp_region.html,
 "temp/sphere"_compute_temp_sphere.html,
 "ti"_compute_ti.html :tb(c=6,ea=c)
 
 These are compute styles contributed by users, which can be used if
 "LAMMPS is built with the appropriate package"_Section_start.html#2_3.
 
 "ackland/atom"_compute_ackland_atom.html,
 "ke/eff"_compute_ke_eff.html,
 "ke/atom/eff"_compute_ke_atom_eff.html,
 "temp/eff"_compute_temp_eff.html,
 "temp/deform/eff"_compute_temp_deform_eff.html,
 "temp/region/eff"_compute_temp_region_eff.html :tb(c=6,ea=c)
 
 :line
 
 Pair_style potentials :h4
 
 See the "pair_style"_pair_style.html command for an overview of pair
 potentials.  Click on the style itself for a full description:
 
 "none"_pair_none.html,
 "hybrid"_pair_hybrid.html,
 "hybrid/overlay"_pair_hybrid.html,
 "airebo"_pair_airebo.html,
 "born"_pair_born.html,
 "born/coul/long"_pair_born.html,
 "buck"_pair_buck.html,
 "buck/coul/cut"_pair_buck.html,
 "buck/coul/long"_pair_buck.html,
 "colloid"_pair_colloid.html,
 "comb"_pair_comb.html,
 "coul/cut"_pair_coul.html,
 "coul/debye"_pair_coul.html,
 "coul/long"_pair_coul.html,
 "dipole/cut"_pair_dipole.html,
 "dpd"_pair_dpd.html,
 "dpd/tstat"_pair_dpd.html,
 "dsmc"_pair_dsmc.html,
 "eam"_pair_eam.html,
 "eam/opt"_pair_eam.html,
 "eam/alloy"_pair_eam.html,
 "eam/alloy/opt"_pair_eam.html,
 "eam/fs"_pair_eam.html,
 "eam/fs/opt"_pair_eam.html,
 "eim"_pair_eim.html,
 "gauss"_pair_gauss.html,
 "gayberne"_pair_gayberne.html,
 "gayberne/gpu"_pair_gayberne.html,
 "gran/hertz/history"_pair_gran.html,
 "gran/hooke"_pair_gran.html,
 "gran/hooke/history"_pair_gran.html,
+"hbond/dreiding/lj"_pair_hbond_dreiding.html,
+"hbond/dreiding/morse"_pair_hbond_dreiding.html,
 "lj/charmm/coul/charmm"_pair_charmm.html,
 "lj/charmm/coul/charmm/implicit"_pair_charmm.html,
 "lj/charmm/coul/long"_pair_charmm.html,
 "lj/charmm/coul/long/opt"_pair_charmm.html,
 "lj/class2"_pair_class2.html,
 "lj/class2/coul/cut"_pair_class2.html,
 "lj/class2/coul/long"_pair_class2.html,
 "lj/cut"_pair_lj.html,
 "lj/cut/gpu"_pair_lj.html,
 "lj/cut/opt"_pair_lj.html,
 "lj/cut/coul/cut"_pair_lj.html,
 "lj/cut/coul/debye"_pair_lj.html,
 "lj/cut/coul/long"_pair_lj.html,
 "lj/cut/coul/long/tip4p"_pair_lj.html,
 "lj/expand"_pair_lj_expand.html,
 "lj/gromacs"_pair_gromacs.html,
 "lj/gromacs/coul/gromacs"_pair_gromacs.html,
 "lj/smooth"_pair_lj_smooth.html,
 "lj96/cut"_pair_lj96_cut.html,
 "lubricate"_pair_lubricate.html,
 "meam"_pair_meam.html,
 "morse"_pair_morse.html,
 "morse/opt"_pair_morse.html,
 "peri/lps"_pair_peri.html,
 "peri/pmb"_pair_peri.html,
 "reax"_pair_reax.html,
 "resquared"_pair_resquared.html,
 "soft"_pair_soft.html,
 "sw"_pair_sw.html,
 "table"_pair_table.html,
 "tersoff"_pair_tersoff.html,
 "tersoff/zbl"_pair_tersoff_zbl.html,
 "yukawa"_pair_yukawa.html,
 "yukawa/colloid"_pair_yukawa_colloid.html :tb(c=4,ea=c)
 
 These are pair styles contributed by users, which can be used if
 "LAMMPS is built with the appropriate package"_Section_start.html#2_3.
 
 "buck/coul"_pair_buck_coul.html,
 "cg/cmm"_pair_cmm.html,
 "cg/cmm/coul/cut"_pair_cmm.html,
 "cg/cmm/coul/long"_pair_cmm.html,
 "eam/cd"_pair_eam.html,
 "eff/cut"_pair_eff.html,
 "lj/coul"_pair_lj_coul.html,
 "reax/c"_pair_reax_c.html :tb(c=4,ea=c)
 
 :line
 
 Bond_style potentials :h4
 
 See the "bond_style"_bond_style.html command for an overview of bond
 potentials.  Click on the style itself for a full description:
 
 "none"_bond_none.html,
 "hybrid"_bond_hybrid.html,
 "class2"_bond_class2.html,
 "fene"_bond_fene.html,
 "fene/expand"_bond_fene_expand.html,
 "harmonic"_bond_harmonic.html,
 "morse"_bond_morse.html,
 "nonlinear"_bond_nonlinear.html,
 "quartic"_bond_quartic.html,
 "table"_bond_table.html :tb(c=4,ea=c,w=100)
 
 :line
 
 Angle_style potentials :h4
 
 See the "angle_style"_angle_style.html command for an overview of
 angle potentials.  Click on the style itself for a full description:
 
 "none"_angle_none.html,
 "hybrid"_angle_hybrid.html,
 "charmm"_angle_charmm.html,
 "class2"_angle_class2.html,
 "cosine"_angle_cosine.html,
 "cosine/delta"_angle_cosine_delta.html,
+"cosine/periodic"_angle_cosine_periodic.html,
 "cosine/squared"_angle_cosine_squared.html,
 "harmonic"_angle_harmonic.html,
 "table"_angle_table.html :tb(c=4,ea=c,w=100)
 
 These are angle styles contributed by users, which can be used if
 "LAMMPS is built with the appropriate package"_Section_start.html#2_3.
 
 "cg/cmm"_angle_cmm.html :tb(c=4,ea=c)
 
 :line
 
 Dihedral_style potentials :h4
 
 See the "dihedral_style"_dihedral_style.html command for an overview
 of dihedral potentials.  Click on the style itself for a full
 description:
 
 "none"_dihedral_none.html,
 "hybrid"_dihedral_hybrid.html,
 "charmm"_dihedral_charmm.html,
 "class2"_dihedral_class2.html,
 "harmonic"_dihedral_harmonic.html,
 "helix"_dihedral_helix.html,
 "multi/harmonic"_dihedral_multi_harmonic.html,
 "opls"_dihedral_opls.html :tb(c=4,ea=c,w=100)
 
 :line
 
 Improper_style potentials :h4
 
 See the "improper_style"_improper_style.html command for an overview
 of improper potentials.  Click on the style itself for a full
 description:
 
 "none"_improper_none.html,
 "hybrid"_improper_hybrid.html,
 "class2"_improper_class2.html,
 "cvff"_improper_cvff.html,
-"harmonic"_improper_harmonic.html :tb(c=4,ea=c,w=100)
+"harmonic"_improper_harmonic.html,
+"umbrella"_improper_umbrella.html :tb(c=4,ea=c,w=100)
 
 :line
 
 Kspace solvers :h4
 
 See the "kspace_style"_kspace_style.html command for an overview of
 Kspace solvers.  Click on the style itself for a full description:
 
 "ewald"_kspace_style.html,
 "pppm"_kspace_style.html,
 "pppm/tip4p"_kspace_style.html :tb(c=4,ea=c,w=100)
 
 These are Kspace solvers contributed by users, which can be used if
 "LAMMPS is built with the appropriate package"_Section_start.html#2_3.
 
 "ewald/n"_kspace_style.html :tb(c=4,ea=c,w=100)
diff --git a/doc/Section_howto.html b/doc/Section_howto.html
index 6fc10741c..0924aab3f 100644
--- a/doc/Section_howto.html
+++ b/doc/Section_howto.html
@@ -1,1729 +1,1691 @@
 <HTML>
 <CENTER><A HREF = "Section_commands.html">Previous Section</A> - <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> - <A HREF = "Section_example.html">Next Section</A> 
 </CENTER>
 
 
 
 
 
 
 <HR>
 
 <H3>4. How-to discussions 
 </H3>
 <P>The following sections describe what commands can be used to perform
 certain kinds of LAMMPS simulations.
 </P>
 4.1 <A HREF = "#4_1">Restarting a simulation</A><BR>
 4.2 <A HREF = "#4_2">2d simulations</A><BR>
-4.3 <A HREF = "#4_3">CHARMM and AMBER force fields</A><BR>
+4.3 <A HREF = "#4_3">CHARMM, AMBER, and DREIDING force fields</A><BR>
 4.4 <A HREF = "#4_4">Running multiple simulations from one input script</A><BR>
 4.5 <A HREF = "#4_5">Multi-replica simulations</A><BR>
 4.6 <A HREF = "#4_6">Granular models</A><BR>
 4.7 <A HREF = "#4_7">TIP3P water model</A><BR>
 4.8 <A HREF = "#4_8">TIP4P water model</A><BR>
 4.9 <A HREF = "#4_9">SPC water model</A><BR>
 4.10 <A HREF = "#4_10">Coupling LAMMPS to other codes</A><BR>
 4.11 <A HREF = "#4_11">Visualizing LAMMPS snapshots</A><BR>
 4.12 <A HREF = "#4_12">Triclinic (non-orthogonal) simulation boxes</A><BR>
 4.13 <A HREF = "#4_13">NEMD simulations</A><BR>
 4.14 <A HREF = "#4_14">Extended spherical and aspherical particles</A><BR>
 4.15 <A HREF = "#4_15">Output from LAMMPS (thermo, dumps, computes, fixes, variables)</A><BR>
 4.16 <A HREF = "#4_16">Thermostatting, barostatting and computing temperature</A><BR>
 4.17 <A HREF = "#4_17">Walls</A><BR>
-4.18 <A HREF = "#4_18">Elastic constants</A><BR>
-4.19 <A HREF = "#4_19">Computing free energies from thermodyanmic integration</A> <BR>
+4.18 <A HREF = "#4_18">Elastic constants</A> <BR>
 
 <P>The example input scripts included in the LAMMPS distribution and
 highlighted in <A HREF = "Section_example.html">this section</A> also show how to
 setup and run various kinds of problems.
 </P>
 <HR>
 
 <A NAME = "4_1"></A><H4>4.1 Restarting a simulation 
 </H4>
 <P>There are 3 ways to continue a long LAMMPS simulation.  Multiple
 <A HREF = "run.html">run</A> commands can be used in the same input script.  Each
 run will continue from where the previous run left off.  Or binary
 restart files can be saved to disk using the <A HREF = "restart.html">restart</A>
 command.  At a later time, these binary files can be read via a
 <A HREF = "read_restart.html">read_restart</A> command in a new script.  Or they can
 be converted to text data files and read by a
 <A HREF = "read_data.html">read_data</A> command in a new script.  <A HREF = "Section_tools.html">This
 section</A> discusses the <I>restart2data</I> tool that is
 used to perform the conversion.
 </P>
 <P>Here we give examples of 2 scripts that read either a binary restart
 file or a converted data file and then issue a new run command to
 continue where the previous run left off.  They illustrate what
 settings must be made in the new script.  Details are discussed in the
 documentation for the <A HREF = "read_restart.html">read_restart</A> and
 <A HREF = "read_data.html">read_data</A> commands.
 </P>
 <P>Look at the <I>in.chain</I> input script provided in the <I>bench</I> directory
 of the LAMMPS distribution to see the original script that these 2
 scripts are based on.  If that script had the line
 </P>
 <PRE>restart	        50 tmp.restart 
 </PRE>
 <P>added to it, it would produce 2 binary restart files (tmp.restart.50
 and tmp.restart.100) as it ran.
 </P>
 <P>This script could be used to read the 1st restart file and re-run the
 last 50 timesteps:
 </P>
 <PRE>read_restart	tmp.restart.50 
 </PRE>
 <PRE>neighbor	0.4 bin
 neigh_modify	every 1 delay 1 
 </PRE>
 <PRE>fix		1 all nve
 fix		2 all langevin 1.0 1.0 10.0 904297 
 </PRE>
 <PRE>timestep	0.012 
 </PRE>
 <PRE>run		50 
 </PRE>
 <P>Note that the following commands do not need to be repeated because
 their settings are included in the restart file: <I>units, atom_style,
 special_bonds, pair_style, bond_style</I>.  However these commands do
 need to be used, since their settings are not in the restart file:
 <I>neighbor, fix, timestep</I>.
 </P>
 <P>If you actually use this script to perform a restarted run, you will
 notice that the thermodynamic data match at step 50 (if you also put a
 "thermo 50" command in the original script), but do not match at step
 100.  This is because the <A HREF = "fix_langevin.html">fix langevin</A> command
 uses random numbers in a way that does not allow for perfect restarts.
 </P>
 <P>As an alternate approach, the restart file could be converted to a data
 file using this tool:
 </P>
 <PRE>restart2data tmp.restart.50 tmp.restart.data 
 </PRE>
 <P>Then, this script could be used to re-run the last 50 steps:
 </P>
 <PRE>units		lj
 atom_style	bond
 pair_style	lj/cut 1.12
 pair_modify	shift yes
 bond_style	fene
 special_bonds   0.0 1.0 1.0 
 </PRE>
 <PRE>read_data	tmp.restart.data 
 </PRE>
 <PRE>neighbor	0.4 bin
 neigh_modify	every 1 delay 1 
 </PRE>
 <PRE>fix		1 all nve
 fix		2 all langevin 1.0 1.0 10.0 904297 
 </PRE>
 <PRE>timestep	0.012 
 </PRE>
 <PRE>reset_timestep	50
 run		50 
 </PRE>
 <P>Note that nearly all the settings specified in the original <I>in.chain</I>
 script must be repeated, except the <I>pair_coeff</I> and <I>bond_coeff</I>
 commands since the new data file lists the force field coefficients.
 Also, the <A HREF = "reset_timestep.html">reset_timestep</A> command is used to tell
 LAMMPS the current timestep.  This value is stored in restart files,
 but not in data files.
 </P>
 <HR>
 
 <A NAME = "4_2"></A><H4>4.2 2d simulations 
 </H4>
 <P>Use the <A HREF = "dimension.html">dimension</A> command to specify a 2d simulation.
 </P>
 <P>Make the simulation box periodic in z via the <A HREF = "boundary.html">boundary</A>
 command.  This is the default.
 </P>
 <P>If using the <A HREF = "create_box.html">create box</A> command to define a
 simulation box, set the z dimensions narrow, but finite, so that the
 create_atoms command will tile the 3d simulation box with a single z
 plane of atoms - e.g.
 </P>
 <PRE><A HREF = "create_box.html">create box</A> 1 -10 10 -10 10 -0.25 0.25 
 </PRE>
 <P>If using the <A HREF = "read_data.html">read data</A> command to read in a file of
 atom coordinates, set the "zlo zhi" values to be finite but narrow,
 similar to the create_box command settings just described.  For each
 atom in the file, assign a z coordinate so it falls inside the
 z-boundaries of the box - e.g. 0.0.
 </P>
 <P>Use the <A HREF = "fix_enforce2d.html">fix enforce2d</A> command as the last
 defined fix to insure that the z-components of velocities and forces
 are zeroed out every timestep.  The reason to make it the last fix is
 so that any forces induced by other fixes will be zeroed out.
 </P>
 <P>Many of the example input scripts included in the LAMMPS distribution
 are for 2d models.
 </P>
 <P>IMPORTANT NOTE: Some models in LAMMPS treat particles as extended
 spheres, as opposed to point particles.  In 2d, the particles will
 still be spheres, not disks, meaning their moment of inertia will be
 the same as in 3d.
 </P>
 <HR>
 
-<A NAME = "4_3"></A><H4>4.3 CHARMM and AMBER force fields 
+<A NAME = "4_3"></A><H4>4.3 CHARMM, AMBER, and DREIDING force fields 
 </H4>
-<P>There are many different ways to compute forces in the <A HREF = "http://www.scripps.edu/brooks">CHARMM</A>
-and <A HREF = "http://amber.scripps.edu">AMBER</A> molecular dynamics codes, only some of which are
-available as options in LAMMPS.  A force field has 2 parts: the
-formulas that define it and the coefficients used for a particular
-system.  Here we only discuss formulas implemented in LAMMPS.  Setting
+<P>A force field has 2 parts: the formulas that define it and the
+coefficients used for a particular system.  Here we only discuss
+formulas implemented in LAMMPS that correspond to formulas commonly
+used in the CHARMM, AMBER, and DREIDING force fields.  Setting
 coefficients is done in the input data file via the
 <A HREF = "read_data.html">read_data</A> command or in the input script with
 commands like <A HREF = "pair_coeff.html">pair_coeff</A> or
-<A HREF = "bond_coeff.html">bond_coeff</A>.  See <A HREF = "Section_tools.html">this section</A> for
-additional tools that can use CHARMM or AMBER to assign force field
-coefficients and convert their output into LAMMPS input.
+<A HREF = "bond_coeff.html">bond_coeff</A>.  See <A HREF = "Section_tools.html">this section</A>
+for additional tools that can use CHARMM or AMBER to assign force
+field coefficients and convert their output into LAMMPS input.
 </P>
 <P>See <A HREF = "#MacKerell">(MacKerell)</A> for a description of the CHARMM force
 field.  See <A HREF = "#Cornell">(Cornell)</A> for a description of the AMBER force
 field.
 </P>
 
 
 
 
 <P>These style choices compute force field formulas that are consistent
 with common options in CHARMM or AMBER.  See each command's
 documentation for the formula it computes.
 </P>
-<UL><LI><A HREF = "bond_style.html">bond_style</A> harmonic
-<LI><A HREF = "angle_style.html">angle_style</A> charmm
-<LI><A HREF = "dihedral_style.html">dihedral_style</A> charmm
-<LI><A HREF = "pair_style.html">pair_style</A> lj/charmm/coul/charmm
-<LI><A HREF = "pair_style.html">pair_style</A> lj/charmm/coul/charmm/implicit
-<LI><A HREF = "pair_style.html">pair_style</A> lj/charmm/coul/long 
+<UL><LI><A HREF = "bond_harmonic.html">bond_style</A> harmonic
+<LI><A HREF = "angle_charmm.html">angle_style</A> charmm
+<LI><A HREF = "dihedral_charmm.html">dihedral_style</A> charmm
+<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/charmm
+<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/charmm/implicit
+<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/long 
 </UL>
 <UL><LI><A HREF = "special_bonds.html">special_bonds</A> charmm
 <LI><A HREF = "special_bonds.html">special_bonds</A> amber 
 </UL>
+<P>DREIDING is a generic force field developed by the <A HREF = "http://www.wag.caltech.edu">Goddard
+group</A> at Caltech and is useful for
+predicting structures and dynamics of organic, biological and
+main-group inorganic molecules. The philosophy in DREIDING is to use
+general force constants and geometry parameters based on simple
+hybridization considerations, rather than individual force constants
+and geometric parameters that depend on the particular combinations of
+atoms involved in the bond, angle, or torsion terms. DREIDING has an
+<A HREF = "pair_hbond_dreiding.html">explicit hydrogen bond term</A> to describe
+interactions involving a hydrogen atom (H___A) on very electronegative
+atoms (N, O, F).
+</P>
+<P>See <A HREF = "#Mayo">(Mayo)</A> for a description of the DREIDING force field
+</P>
+<P>These style choices compute force field formulas that are consistent
+with the DREIDING force field.  See each command's
+documentation for the formula it computes.
+</P>
+<UL><LI><A HREF = "bond_harmonic.html">bond_style</A> harmonic
+<LI><A HREF = "bond_morse.html">bond_style</A> morse 
+</UL>
+<UL><LI><A HREF = "angle_harmonic.html">angle_style</A> harmonic
+<LI><A HREF = "angle_cosine.html">angle_style</A> cosine
+<LI><A HREF = "angle_cosine_periodic.html">angle_style</A> cosine/periodic 
+</UL>
+<UL><LI><A HREF = "dihedral_charmm.html">dihedral_style</A> charmm
+<LI><A HREF = "improper_umbrella.html">improper_style</A> umbrella 
+</UL>
+<UL><LI><A HREF = "pair_buck.html">pair_style</A> buck
+<LI><A HREF = "pair_buck.html">pair_style</A> buck/coul/cut
+<LI><A HREF = "pair_buck.html">pair_style</A> buck/coul/long
+<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut
+<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut/coul/cut
+<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut/coul/long 
+</UL>
+<UL><LI><A HREF = "pair_hbond_dreiding.html">pair_style</A> hbond/dreiding/lj
+<LI><A HREF = "pair_hbond_dreiding.html">pair_style</A> hbond/dreiding/morse 
+</UL>
+<UL><LI><A HREF = "special_bonds.html">special_bonds</A> dreiding 
+</UL>
 <HR>
 
 <A NAME = "4_4"></A><H4>4.4 Running multiple simulations from one input script 
 </H4>
 <P>This can be done in several ways.  See the documentation for
 individual commands for more details on how these examples work.
 </P>
 <P>If "multiple simulations" means continue a previous simulation for
 more timesteps, then you simply use the <A HREF = "run.html">run</A> command
 multiple times.  For example, this script
 </P>
 <PRE>units lj
 atom_style atomic
 read_data data.lj
 run 10000
 run 10000
 run 10000
 run 10000
 run 10000 
 </PRE>
 <P>would run 5 successive simulations of the same system for a total of
 50,000 timesteps.
 </P>
 <P>If you wish to run totally different simulations, one after the other,
 the <A HREF = "clear.html">clear</A> command can be used in between them to
 re-initialize LAMMPS.  For example, this script
 </P>
 <PRE>units lj
 atom_style atomic
 read_data data.lj
 run 10000
 clear
 units lj
 atom_style atomic
 read_data data.lj.new
 run 10000 
 </PRE>
 <P>would run 2 independent simulations, one after the other.
 </P>
 <P>For large numbers of independent simulations, you can use
 <A HREF = "variable.html">variables</A> and the <A HREF = "next.html">next</A> and
 <A HREF = "jump.html">jump</A> commands to loop over the same input script
 multiple times with different settings.  For example, this
 script, named in.polymer
 </P>
 <PRE>variable d index run1 run2 run3 run4 run5 run6 run7 run8
 shell cd $d
 read_data data.polymer
 run 10000
 shell cd ..
 clear
 next d
 jump in.polymer 
 </PRE>
 <P>would run 8 simulations in different directories, using a data.polymer
 file in each directory.  The same concept could be used to run the
 same system at 8 different temperatures, using a temperature variable
 and storing the output in different log and dump files, for example
 </P>
 <PRE>variable a loop 8
 variable t index 0.8 0.85 0.9 0.95 1.0 1.05 1.1 1.15
 log log.$a
 read data.polymer
 velocity all create $t 352839
 fix 1 all nvt $t $t 100.0
 dump 1 all atom 1000 dump.$a
 run 100000
 next t
 next a
 jump in.polymer 
 </PRE>
 <P>All of the above examples work whether you are running on 1 or
 multiple processors, but assumed you are running LAMMPS on a single
 partition of processors.  LAMMPS can be run on multiple partitions via
 the "-partition" command-line switch as described in <A HREF = "Section_start.html#2_6">this
 section</A> of the manual.
 </P>
 <P>In the last 2 examples, if LAMMPS were run on 3 partitions, the same
 scripts could be used if the "index" and "loop" variables were
 replaced with <I>universe</I>-style variables, as described in the
 <A HREF = "variable.html">variable</A> command.  Also, the "next t" and "next a"
 commands would need to be replaced with a single "next a t" command.
 With these modifications, the 8 simulations of each script would run
 on the 3 partitions one after the other until all were finished.
 Initially, 3 simulations would be started simultaneously, one on each
 partition.  When one finished, that partition would then start
 the 4th simulation, and so forth, until all 8 were completed.
 </P>
 <HR>
 
 <A NAME = "4_5"></A><H4>4.5 Multi-replica simulations 
 </H4>
 <P>Several commands in LAMMPS run mutli-replica simulations, meaning
 that multiple instances (replicas) of your simulation are run
 simultaneously, with small amounts of data exchanged between replicas
 periodically.
 </P>
 <P>These are the relevant commands:
 </P>
 <UL><LI><A HREF = "neb.html">neb</A> for nudged elastic band calculations
 <LI><A HREF = "prd.html">prd</A> for parallel replica dynamics
 <LI><A HREF = "temper.html">temper</A> for parallel tempering 
 </UL>
 <P>NEB is a method for finding transition states and barrier energies.
 PRD is a method for performing accelerated dynamics to find and
 perform infrequent events.  Parallel tempering or replica exchange
 runs different replicas at a series of temperature to facilitate
 rare-event sampling.
 </P>
 <P>In all these cases, you must run with one or more processors per
 replica.  The processors assigned to each replica are determined at
 run-time by using the <A HREF = "Section_start.html#2_6">-partition command-line
 switch</A> to launch LAMMPS on multiple
 partitions, which in this context are the same as replicas.  E.g.
 these commands:
 </P>
 <PRE>mpirun -np 16 lmp_linux -partition 8x2 -in in.temper
 mpirun -np 8 lmp_linux -partition 8x1 -in in.neb 
 </PRE>
 <P>would each run 8 replicas, on either 16 or 8 processors.  Note the use
 of the <A HREF = "Section_start.html#2_6">-in command-line switch</A> to specify the
 input script which is required when running in multi-replica mode.
 </P>
 <P>Also note that with MPI installed on a machine (e.g. your desktop),
 you can run on more (virtual) processors than you have physical
 processors.  Thus the above commands could be run on a
 single-processor (or few-processor) desktop so that you can run
 a multi-replica simulation on more replicas than you have
 physical processors.
 </P>
 <HR>
 
 <A NAME = "4_6"></A><H4>4.6 Granular models 
 </H4>
 <P>Granular system are composed of spherical particles with a diameter,
 as opposed to point particles.  This means they have an angular
 velocity and torque can be imparted to them to cause them to rotate.
 </P>
 <P>To run a simulation of a granular model, you will want to use
 the following commands:
 </P>
 <UL><LI><A HREF = "atom_style.html">atom_style</A> granular
 <LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
 <LI><A HREF = "fix_gravity.html">fix gravity</A> 
 </UL>
 <P>This compute
 </P>
 <UL><LI><A HREF = "compute_erotate_sphere.html">compute erotate/sphere</A> 
 </UL>
 <P>calculates rotational kinetic energy which can be <A HREF = "Section_howto.html#4_15">output with
 thermodynamic info</A>.
 </P>
 <P>Use one of these 3 pair potentials, which compute forces and torques
 between interacting pairs of particles:
 </P>
 <UL><LI><A HREF = "pair_style.html">pair_style</A> gran/history
 <LI><A HREF = "pair_style.html">pair_style</A> gran/no_history
 <LI><A HREF = "pair_style.html">pair_style</A> gran/hertzian 
 </UL>
 <P>These commands implement fix options specific to granular systems:
 </P>
 <UL><LI><A HREF = "fix_freeze.html">fix freeze</A>
 <LI><A HREF = "fix_pour.html">fix pour</A>
 <LI><A HREF = "fix_viscous.html">fix viscous</A>
 <LI><A HREF = "fix_wall_gran.html">fix wall/gran</A> 
 </UL>
 <P>The fix style <I>freeze</I> zeroes both the force and torque of frozen
 atoms, and should be used for granular system instead of the fix style
 <I>setforce</I>.
 </P>
 <P>For computational efficiency, you can eliminate needless pairwise
 computations between frozen atoms by using this command:
 </P>
 <UL><LI><A HREF = "neigh_modify.html">neigh_modify</A> exclude 
 </UL>
 <HR>
 
 <A NAME = "4_7"></A><H4>4.7 TIP3P water model 
 </H4>
 <P>The TIP3P water model as implemented in CHARMM
 <A HREF = "#MacKerell">(MacKerell)</A> specifies a 3-site rigid water molecule with
 charges and Lennard-Jones parameters assigned to each of the 3 atoms.
 In LAMMPS the <A HREF = "fix_shake.html">fix shake</A> command can be used to hold
 the two O-H bonds and the H-O-H angle rigid.  A bond style of
 <I>harmonic</I> and an angle style of <I>harmonic</I> or <I>charmm</I> should also be
 used.
 </P>
 <P>These are the additional parameters (in real units) to set for O and H
 atoms and the water molecule to run a rigid TIP3P-CHARMM model with a
 cutoff.  The K values can be used if a flexible TIP3P model (without
 fix shake) is desired.  If the LJ epsilon and sigma for HH and OH are
 set to 0.0, it corresponds to the original 1983 TIP3P model
 <A HREF = "#Jorgensen">(Jorgensen)</A>.
 </P>
 <P>O mass = 15.9994<BR>
 H mass = 1.008 <BR>
 </P>
 <P>O charge = -0.834<BR>
 H charge = 0.417 <BR>
 </P>
 <P>LJ epsilon of OO = 0.1521<BR>
 LJ sigma of OO = 3.1507<BR>
 LJ epsilon of HH = 0.0460<BR>
 LJ sigma of HH = 0.4000<BR>
 LJ epsilon of OH = 0.0836<BR>
 LJ sigma of OH = 1.7753 <BR>
 </P>
 <P>K of OH bond = 450<BR>
 r0 of OH bond = 0.9572 <BR>
 </P>
 <P>K of HOH angle = 55<BR>
 theta of HOH angle = 104.52 <BR>
 </P>
 <P>These are the parameters to use for TIP3P with a long-range Coulombic
 solver (Ewald or PPPM in LAMMPS), see <A HREF = "#Price">(Price)</A> for details:
 </P>
 <P>O mass = 15.9994<BR>
 H mass = 1.008 <BR>
 </P>
 <P>O charge = -0.830<BR>
 H charge = 0.415 <BR>
 </P>
 <P>LJ epsilon of OO = 0.102<BR>
 LJ sigma of OO = 3.188<BR>
 LJ epsilon, sigma of OH, HH = 0.0 <BR>
 </P>
 <P>K of OH bond = 450<BR>
 r0 of OH bond = 0.9572 <BR>
 </P>
 <P>K of HOH angle = 55<BR>
 theta of HOH angle = 104.52 <BR>
 </P>
 <P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
 models</A>.
 </P>
 <HR>
 
 <A NAME = "4_8"></A><H4>4.8 TIP4P water model 
 </H4>
 <P>The four-point TIP4P rigid water model extends the traditional
 three-point TIP3P model by adding an additional site, usually
 massless, where the charge associated with the oxygen atom is placed.
 This site M is located at a fixed distance away from the oxygen along
 the bisector of the HOH bond angle.  A bond style of <I>harmonic</I> and an
 angle style of <I>harmonic</I> or <I>charmm</I> should also be used.
 </P>
 <P>Currently, only a four-point model for long-range Coulombics is
 implemented via the LAMMPS <A HREF = "pair_lj.html">pair style
 lj/cut/coul/long/tip4p</A>.  A cutoff version may be added
 the future.  For both models, the bond lengths and bond angles should
 be held fixed using the <A HREF = "fix_shake.html">fix shake</A> command.
 </P>
 <P>These are the additional parameters (in real units) to set for O and H
 atoms and the water molecule to run a rigid TIP4P model with a cutoff
 <A HREF = "#Jorgensen">(Jorgensen)</A>.  Note that the OM distance is specified in
 the <A HREF = "pair_style.html">pair_style</A> command, not as part of the pair
 coefficients.
 </P>
 <P>O mass = 15.9994<BR>
 H mass = 1.008 <BR>
 </P>
 <P>O charge = -1.040<BR>
 H charge = 0.520 <BR>
 </P>
 <P>r0 of OH bond = 0.9572<BR>
 theta of HOH angle = 104.52 <BR>
 </P>
 <P>OM distance = 0.15 <BR>
 </P>
 <P>LJ epsilon of O-O = 0.1550<BR>
 LJ sigma of O-O = 3.1536<BR>
 LJ epsilon, sigma of OH, HH = 0.0 <BR>
 </P>
 <P>These are the parameters to use for TIP4P with a long-range Coulombic
 solver (Ewald or PPPM in LAMMPS):
 </P>
 <P>O mass = 15.9994<BR>
 H mass = 1.008 <BR>
 </P>
 <P>O charge = -1.0484<BR>
 H charge = 0.5242 <BR>
 </P>
 <P>r0 of OH bond = 0.9572<BR>
 theta of HOH angle = 104.52 <BR>
 </P>
 <P>OM distance = 0.1250 <BR>
 </P>
 <P>LJ epsilon of O-O = 0.16275<BR>
 LJ sigma of O-O = 3.16435<BR>
 LJ epsilon, sigma of OH, HH = 0.0 <BR>
 </P>
 <P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
 models</A>.
 </P>
 <HR>
 
 <A NAME = "4_9"></A><H4>4.9 SPC water model 
 </H4>
 <P>The SPC water model specifies a 3-site rigid water molecule with
 charges and Lennard-Jones parameters assigned to each of the 3 atoms.
 In LAMMPS the <A HREF = "fix_shake.html">fix shake</A> command can be used to hold
 the two O-H bonds and the H-O-H angle rigid.  A bond style of
 <I>harmonic</I> and an angle style of <I>harmonic</I> or <I>charmm</I> should also be
 used.
 </P>
 <P>These are the additional parameters (in real units) to set for O and H
 atoms and the water molecule to run a rigid SPC model.
 </P>
 <P>O mass = 15.9994<BR>
 H mass = 1.008 <BR>
 </P>
 <P>O charge = -0.820<BR>
 H charge = 0.410 <BR>
 </P>
 <P>LJ epsilon of OO = 0.1553<BR>
 LJ sigma of OO = 3.166<BR>
 LJ epsilon, sigma of OH, HH = 0.0 <BR>
 </P>
 <P>r0 of OH bond = 1.0<BR>
 theta of HOH angle = 109.47 <BR>
 </P>
 <P>Note that as originally proposed, the SPC model was run with a 9
 Angstrom cutoff for both LJ and Coulommbic terms.  It can also be used
 with long-range Coulombics (Ewald or PPPM in LAMMPS), without changing
 any of the parameters above, though it becomes a different model in
 that mode of usage.
 </P>
 <P>The SPC/E (extended) water model is the same, except
 the partial charge assignemnts change:
 </P>
 <P>O charge = -0.8476<BR>
 H charge = 0.4238 <BR>
 </P>
 <P>See the <A HREF = "#Berendsen">(Berendsen)</A> reference for more details on both
 the SPC and SPC/E models.
 </P>
 <P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
 models</A>.
 </P>
 <HR>
 
 <A NAME = "4_10"></A><H4>4.10 Coupling LAMMPS to other codes 
 </H4>
 <P>LAMMPS is designed to allow it to be coupled to other codes.  For
 example, a quantum mechanics code might compute forces on a subset of
 atoms and pass those forces to LAMMPS.  Or a continuum finite element
 (FE) simulation might use atom positions as boundary conditions on FE
 nodal points, compute a FE solution, and return interpolated forces on
 MD atoms.
 </P>
 <P>LAMMPS can be coupled to other codes in at least 3 ways.  Each has
 advantages and disadvantages, which you'll have to think about in the
 context of your application.
 </P>
 <P>(1) Define a new <A HREF = "fix.html">fix</A> command that calls the other code.  In
 this scenario, LAMMPS is the driver code.  During its timestepping,
 the fix is invoked, and can make library calls to the other code,
 which has been linked to LAMMPS as a library.  This is the way the
 <A HREF = "http://www.rpi.edu/~anderk5/lab">POEMS</A> package that performs constrained rigid-body motion on
 groups of atoms is hooked to LAMMPS.  See the
 <A HREF = "fix_poems.html">fix_poems</A> command for more details.  See <A HREF = "Section_modify.html">this
 section</A> of the documentation for info on how to add
 a new fix to LAMMPS.
 </P>
 
 
 <P>(2) Define a new LAMMPS command that calls the other code.  This is
 conceptually similar to method (1), but in this case LAMMPS and the
 other code are on a more equal footing.  Note that now the other code
 is not called during the timestepping of a LAMMPS run, but between
 runs.  The LAMMPS input script can be used to alternate LAMMPS runs
 with calls to the other code, invoked via the new command.  The
 <A HREF = "run.html">run</A> command facilitates this with its <I>every</I> option, which
 makes it easy to run a few steps, invoke the command, run a few steps,
 invoke the command, etc.
 </P>
 <P>In this scenario, the other code can be called as a library, as in
 (1), or it could be a stand-alone code, invoked by a system() call
 made by the command (assuming your parallel machine allows one or more
 processors to start up another program).  In the latter case the
 stand-alone code could communicate with LAMMPS thru files that the
 command writes and reads.
 </P>
 <P>See <A HREF = "Section_modify.html">this section</A> of the documentation for how to
 add a new command to LAMMPS.
 </P>
 <P>(3) Use LAMMPS as a library called by another code.  In this case the
 other code is the driver and calls LAMMPS as needed.  Or a wrapper
 code could link and call both LAMMPS and another code as libraries.
 Again, the <A HREF = "run.html">run</A> command has options that allow it to be
 invoked with minimal overhead (no setup or clean-up) if you wish to do
 multiple short runs, driven by another program.
 </P>
 <P>Examples of driver codes that call LAMMPS as a library are included in
 the "couple" directory of the LAMMPS distribution; see couple/README
 for more details:
 </P>
 <UL><LI>simple: simple driver programs in C++ and C which invoke LAMMPS as a
 library 
 
 <LI>lammps_quest: coupling of LAMMPS and <A HREF = "http://dft.sandia.gov/Quest">Quest</A>, to run classical
 MD with quantum forces calculated by a density functional code 
 
 <LI>lammps_spparks: coupling of LAMMPS and <A HREF = "http://www.sandia.gov/~sjplimp/spparks.html">SPPARKS</A>, to couple
 a kinetic Monte Carlo model for grain growth using MD to calculate
 strain induced across grain boundaries 
 </UL>
 
 
 
 
 <P><A HREF = "Section_start.html#2_4">This section</A> of the documentation describes
 how to build LAMMPS as a library.  Once this is done, you can
 interface with LAMMPS either via C++, C, or Fortran (or any other
 language that supports a vanilla C-like interface, e.g. a scripting
 language).  For example, from C++ you could create one (or more)
 "instances" of LAMMPS, pass it an input script to process, or execute
 individual commands, all by invoking the correct class methods in
 LAMMPS.  From C or Fortran you can make function calls to do the same
 things.  Library.cpp and library.h contain such a C interface with the
 functions:
 </P>
 <PRE>void lammps_open(int, char **, MPI_Comm, void **);
 void lammps_close(void *);
 void lammps_file(void *, char *);
 char *lammps_command(void *, char *); 
 </PRE>
 <P>The functions contain C++ code you could write in a C++ application
 that was invoking LAMMPS directly.  Note that LAMMPS classes are
 defined within a LAMMPS namespace (LAMMPS_NS) if you use them
 from another C++ application.
 </P>
 <P>Two of the routines in library.cpp are of particular note.  The
 lammps_open() function initiates LAMMPS and takes an MPI communicator
 as an argument.  It returns a pointer to a LAMMPS "object".  As with
 C++, the lammps_open() function can be called multiple times, to
 create multiple instances of LAMMPS.
 </P>
 <P>LAMMPS will run on the set of processors in the communicator.  This
 means the calling code can run LAMMPS on all or a subset of
 processors.  For example, a wrapper script might decide to alternate
 between LAMMPS and another code, allowing them both to run on all the
 processors.  Or it might allocate half the processors to LAMMPS and
 half to the other code and run both codes simultaneously before
 syncing them up periodically.
 </P>
 <P>Library.cpp contains a lammps_command() function to which the caller
 passes a single LAMMPS command (a string).  Thus the calling code can
 read or generate a series of LAMMPS commands (e.g. an input script)
 one line at a time and pass it thru the library interface to setup a
 problem and then run it.
 </P>
 <P>A few other sample functions are included in library.cpp, but the key
 idea is that you can write any functions you wish to define an
 interface for how your code talks to LAMMPS and add them to
 library.cpp and library.h.  The routines you add can access any LAMMPS
 data.  The examples/couple directory has example C++ and C codes which
 show how a stand-alone code can link LAMMPS as a library, run LAMMPS
 on a subset of processors, grab data from LAMMPS, change it, and put
 it back into LAMMPS.
 </P>
 <HR>
 
 <A NAME = "4_11"></A><H4>4.11 Visualizing LAMMPS snapshots 
 </H4>
 <P>LAMMPS itself does not do visualization, but snapshots from LAMMPS
 simulations can be visualized (and analyzed) in a variety of ways.
 </P>
 <P>LAMMPS snapshots are created by the <A HREF = "dump.html">dump</A> command which can
 create files in several formats.  The native LAMMPS dump format is a
 text file (see "dump atom" or "dump custom") which can be visualized
 by the <A HREF = "Section_tools.html#xmovie">xmovie</A> program, included with the
 LAMMPS package.  This produces simple, fast 2d projections of 3d
 systems, and can be useful for rapid debugging of simulation geometry
 and atom trajectories.
 </P>
 <P>Several programs included with LAMMPS as auxiliary tools can convert
 native LAMMPS dump files to other formats.  See the
 <A HREF = "Section_tools.html">Section_tools</A> doc page for details.  The first is
 the <A HREF = "Section_tools.html#charmm">ch2lmp tool</A>, which contains a
 lammps2pdb Perl script which converts LAMMPS dump files into PDB
 files.  The second is the <A HREF = "Section_tools.html#arc">lmp2arc tool</A> which
 converts LAMMPS dump files into Accelrys' Insight MD program files.
 The third is the <A HREF = "Section_tools.html#cfg">lmp2cfg tool</A> which converts
 LAMMPS dump files into CFG files which can be read into the
 <A HREF = "http://mt.seas.upenn.edu/Archive/Graphics/A">AtomEye</A> visualizer.
 </P>
 <P>A Python-based toolkit distributed by our group can read native LAMMPS
 dump files, including custom dump files with additional columns of
 user-specified atom information, and convert them to various formats
 or pipe them into visualization software directly.  See the <A HREF = "http://www.sandia.gov/~sjplimp/pizza.html">Pizza.py
 WWW site</A> for details.  Specifically, Pizza.py can convert
 LAMMPS dump files into PDB, XYZ, <A HREF = "http://www.ensight.com">Ensight</A>, and VTK formats.
 Pizza.py can pipe LAMMPS dump files directly into the Raster3d and
 RasMol visualization programs.  Pizza.py has tools that do interactive
 3d OpenGL visualization and one that creates SVG images of dump file
 snapshots.
 </P>
 <P>LAMMPS can create XYZ files directly (via "dump xyz") which is a
 simple text-based file format used by many visualization programs
 including <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A>.
 </P>
 <P>LAMMPS can create DCD files directly (via "dump dcd") which can be
 read by <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A> in conjunction with a CHARMM PSF file.  Using this
 form of output avoids the need to convert LAMMPS snapshots to PDB
 files.  See the <A HREF = "dump.html">dump</A> command for more information on DCD
 files.
 </P>
 <P>LAMMPS can create XTC files directly (via "dump xtc") which is GROMACS
 file format which can also be read by <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A> for visualization.
 See the <A HREF = "dump.html">dump</A> command for more information on XTC files.
 </P>
 
 
 
 
 
 
 
 
 <HR>
 
 <A NAME = "4_12"></A><H4>4.12 Triclinic (non-orthogonal) simulation boxes 
 </H4>
 <P>By default, LAMMPS uses an orthogonal simulation box to encompass the
 particles.  The <A HREF = "boundary.html">boundary</A> command sets the boundary
 conditions of the box (periodic, non-periodic, etc).  The orthogonal
 box has its "origin" at (xlo,ylo,zlo) and is defined by 3 edge vectors
 starting from the origin given by A = (xhi-xlo,0,0); B =
 (0,yhi-ylo,0); C = (0,0,zhi-zlo).  The 6 parameters
 (xlo,xhi,ylo,yhi,zlo,zhi) are defined at the time the simluation box
 is created, e.g. by the <A HREF = "create_box.html">create_box</A> or
 <A HREF = "read_data.html">read_data</A> or <A HREF = "read_restart.html">read_restart</A>
 commands.  Additionally, LAMMPS defines box size parameters lx,ly,lz
 where lx = xhi-xlo, and similarly in the y and z dimensions.  The 6
 parameters, as well as lx,ly,lz, can be output via the <A HREF = "thermo_style.html">thermo_style
 custom</A> command.
 </P>
 <P>LAMMPS also allows simulations to be perfored in non-orthogonal
 simulation boxes shaped as a parallelepiped with triclinic symmetry.
 The parallelepiped has its "origin" at (xlo,ylo,zlo) and is defined by
 3 edge vectors starting from the origin given by A = (xhi-xlo,0,0); B
 = (xy,yhi-ylo,0); C = (xz,yz,zhi-zlo).  <I>Xy,xz,yz</I> can be 0.0 or
 positive or negative values and are called "tilt factors" because they
 are the amount of displacement applied to faces of an originally
 orthogonal box to transform it into the parallelepiped.  Note that in
 LAMMPS the triclinic simulation box edge vectors A,B,C cannot be
 arbitrary vectors.  As indicated, A must be aligned with the x axis, B
 must be in the xy plane, and C is arbitrary.  However, this is not a
 restriction since it is possible to rotate any set of 3 crystal basis
 vectors so that they meet this restriction.
 </P>
 <P>The 9 parameters (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) are defined at the
 time the simluation box is created.  This happens in one of 3 ways.
 If the <A HREF = "create_box.html">create_box</A> command is used with a region of
 style <I>prism</I>, then a triclinic box is setup.  See the
 <A HREF = "region.html">region</A> command for details.  If the
 <A HREF = "read_data.html">read_data</A> command is used to define the simulation
 box, and the header of the data file contains a line with the "xy xz
 yz" keyword, then a triclinic box is setup.  See the
 <A HREF = "read_data.html">read_data</A> command for details.  Finally, if the
 <A HREF = "read_restart.html">read_restart</A> command reads a restart file which
 was written from a simulation using a triclinic box, then a triclinic
 box will be setup for the restarted simulation.
 </P>
 <P>Note that you can define a triclinic box with all 3 tilt factors =
 0.0, so that it is initially orthogonal.  This is necessary if the box
 will become non-orthogonal, e.g. due to the <A HREF = "fix_nh.html">fix npt</A> or
 <A HREF = "fix_deform.html">fix deform</A> commands.  Alternatively, you can use the
 <A HREF = "change_box.html">change_box</A> command to convert a simulation box from
 orthogonal to triclinic and vice versa.
 </P>
 <P>As with orthogonal boxes, LAMMPS defines triclinic box size parameters
 lx,ly,lz where lx = xhi-xlo, and similarly in the y and z dimensions.
 The 9 parameters, as well as lx,ly,lz, can be output via the
 <A HREF = "thermo_style.html">thermo_style custom</A> command.
 </P>
 <P>To avoid extremely tilted boxes (which would be computationally
 inefficient), no tilt factor can skew the box more than half the
 distance of the parallel box length, which is the 1st dimension in the
 tilt factor (x for xz).  For example, if xlo = 2 and xhi = 12, then
 the x box length is 10 and the xy tilt factor must be between -5 and
 5.  Similarly, both xz and yz must be between -(xhi-xlo)/2 and
 +(yhi-ylo)/2.  Note that this is not a limitation, since if the
 maximum tilt factor is 5 (as in this example), then configurations
 with tilt = ..., -15, -5, 5, 15, 25, ... are geometrically all
 equivalent.
 </P>
 <P>Triclinic crystal structures are often defined using three lattice
 constants <I>a</I>, <I>b</I>, and <I>c</I>, and three angles <I>alpha</I>, <I>beta</I> and
 <I>gamma</I>.  Note that in this nomenclature, the a,b,c lattice constants
 are the scalar lengths of the 3 A,B,C edge vectors defined above.  The
 relationship between these 6 quantities (a,b,c,alpha,beta,gamma) and
 the LAMMPS box sizes (lx,ly,lz) = (xhi-xlo,yhi-ylo,zhi-zlo) and tilt
 factors (xy,xz,yz) is as follows:
 </P>
 <CENTER><IMG SRC = "Eqs/box.jpg">
 </CENTER>
 <P>As discussed on the <A HREF = "dump.html">dump</A> command doc page, when the BOX
 BOUNDS for a snapshot is written to a dump file for a triclinic box,
 an orthogonal bounding box which encloses the triclinic simulation box
 is output, along with the 3 tilt factors (xy, xz, yz) of the triclinic
 box, formatted as follows:
 </P>
 <PRE>ITEM: BOX BOUNDS xy xz yz
 xlo_bound xhi_bound xy
 ylo_bound yhi_bound xz
 zlo_bound zhi_bound yz 
 </PRE>
 <P>This bounding box is convenient for many visualization programs and is
 calculated from the 9 triclinic box parameters
 (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) as follows:
 </P>
 <PRE>xlo_bound = xlo + MIN(0.0,xy,xz,xy+xz)
 xhi_bound = xhi + MAX(0.0,xy,xz,xy+xz)
 ylo_bound = ylo + MIN(0.0,yz)
 yhi_bound = yhi + MAX(0.0,yz)
 zlo_bound = zlo
 zhi_bound = zhi 
 </PRE>
 <P>These formulas can be inverted if you need to convert the bounding box
 back into the triclinic box parameters, e.g. xlo = xlo_bound -
 MIN(0.0,xy,xz,xy+xz).
 </P>
 <P>One use of triclinic simulation boxes is to model solid-state crystals
 with triclinic symmetry.  The <A HREF = "lattice.html">lattice</A> command can be
 used with non-orthogonal basis vectors to define a lattice that will
 tile a triclinic simulation box via the
 <A HREF = "create_atoms.html">create_atoms</A> command.
 </P>
 <P>A second use is to run Parinello-Rahman dyanamics via the <A HREF = "fix_nh.html">fix
 npt</A> command, which will adjust the xy, xz, yz tilt
 factors to compensate for off-diagonal components of the pressure
 tensor.  The analalog for an <A HREF = "minimize.html">energy minimization</A> is
 the <A HREF = "fix_box_relax.html">fix box/relax</A> command.
 </P>
 <P>A third use is to shear a bulk solid to study the response of the
 material.  The <A HREF = "fix_deform.html">fix deform</A> command can be used for
 this purpose.  It allows dynamic control of the xy, xz, yz tilt
 factors as a simulation runs.  This is discussed in the next section
 on non-equilibrium MD (NEMD) simulations.
 </P>
 <HR>
 
 <A NAME = "4_13"></A><H4>4.13 NEMD simulations 
 </H4>
 <P>Non-equilibrium molecular dynamics or NEMD simulations are typically
 used to measure a fluid's rheological properties such as viscosity.
 In LAMMPS, such simulations can be performed by first setting up a
 non-orthogonal simulation box (see the preceding Howto section).
 </P>
 <P>A shear strain can be applied to the simulation box at a desired
 strain rate by using the <A HREF = "fix_deform.html">fix deform</A> command.  The
 <A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> command can be used to thermostat
 the sheared fluid and integrate the SLLOD equations of motion for the
 system.  Fix nvt/sllod uses <A HREF = "compute_temp_deform.html">compute
 temp/deform</A> to compute a thermal temperature
 by subtracting out the streaming velocity of the shearing atoms.  The
 velocity profile or other properties of the fluid can be monitored via
 the <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command.
 </P>
 <P>As discussed in the previous section on non-orthogonal simulation
 boxes, the amount of tilt or skew that can be applied is limited by
 LAMMPS for computational efficiency to be 1/2 of the parallel box
 length.  However, <A HREF = "fix_deform.html">fix deform</A> can continuously strain
 a box by an arbitrary amount.  As discussed in the <A HREF = "fix_deform.html">fix
 deform</A> command, when the tilt value reaches a limit,
 the box is re-shaped to the opposite limit which is an equivalent
 tiling of periodic space.  The strain rate can then continue to change
 as before.  In a long NEMD simulation these box re-shaping events may
 occur many times.
 </P>
 <P>In a NEMD simulation, the "remap" option of <A HREF = "fix_deform.html">fix
 deform</A> should be set to "remap v", since that is what
 <A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> assumes to generate a velocity
 profile consistent with the applied shear strain rate.
 </P>
 <P>An alternative method for calculating viscosities is provided via the
 <A HREF = "fix_viscosity.html">fix viscosity</A> command.
 </P>
 <HR>
 
 <A NAME = "4_14"></A><H4>4.14 Extended spherical and aspherical particles 
 </H4>
 <P>Typical MD models treat atoms or particles as point masses.
 Sometimes, however, it is desirable to have a model with finite-size
 particles such as spherioids or aspherical ellipsoids.  The difference
 is that such particles have a moment of inertia, rotational energy,
 and angular momentum.  Rotation is induced by torque from interactions
 with other particles.
 </P>
 <P>LAMMPS has several options for running simulations with these kinds of
 particles.  The following aspects are discussed in turn:
 </P>
 <UL><LI>atom styles
 <LI>pair potentials
 <LI>time integration
 <LI>computes, thermodynamics, and dump output
 <LI>rigid bodies composed of extended particles 
 </UL>
 <H5>Atom styles 
 </H5>
 <P>There are 3 <A HREF = "atom_style.html">atom styles</A> that allow for definition of
 finite-size particles: granular, dipole, ellipsoid.
 </P>
 <P>Granular particles are spheriods and each particle can have a unique
 diameter and mass (or density).  These particles store an angular
 velocity (omega) and can be acted upon by torque.
 </P>
 <P>Dipolar particles are typically spheriods with a point dipole and each
 particle type has a diamater and mass, set by the <A HREF = "shape.html">shape</A>
 and <A HREF = "mass.html">mass</A> commands.  These particles store an angular
 velocity (omega) and can be acted upon by torque.  They also store an
 orientation for the point dipole (mu) which has a length set by the
 <A HREF = "dipole.html">dipole</A> command.  The <A HREF = "set.html">set</A> command can be used
 to initialize the orientation of dipole moments.
 </P>
 <P>Ellipsoid particles are aspherical.  Each particle type has an
 ellipsoidal shape and mass, defined by the <A HREF = "shape.html">shape</A> and
 <A HREF = "mass.html">mass</A> commands.  These particles store an angular momentum
 and their orientation (quaternion), and can be acted upon by torque.
 They do not store an angular velocity (omega), which can be in a
 different direction than angular momentum, rather they compute it as
 needed.  Ellipsoidal particles can also store a dipole moment if an
 <A HREF = "atom_style.html">atom_style hybrid ellipsoid dipole</A> is used.  The
 <A HREF = "set.html">set</A> command can be used to initialize the orientation of
 ellipsoidal particles and has a brief explanation of quaternions.
 </P>
 <P>Note that if one of these atom styles is used (or multiple styles via
 the <A HREF = "atom_style.html">atom_style hybrid</A> command), not all particles in
 the system are required to be finite-size or aspherical.  For example,
 if the 3 shape parameters are set to the same value, the particle will
 be a spheroid rather than an ellipsoid.  If the 3 shape parameters are
 all set to 0.0 or if the diameter is set to 0.0, it will be a point
 particle.  If the dipole moment is set to zero, the particle will not
 have a point dipole associated with it.  The pair styles used to
 compute pairwise interactions will typically compute the correct
 interaction in these simplified (cheaper) cases.  <A HREF = "pair_hybrid.html">Pair_style
 hybrid</A> can be used to insure the correct
 interactions are computed for the appropriate style of interactions.
 Likewise, using groups to partition particles (ellipsoid versus
 spheroid versus point particles) will allow you to use the appropriate
 time integrators and temperature computations for each class of
 particles.  See the doc pages for various commands for details.
 </P>
 <P>Also note that for <A HREF = "dimension.html">2d simulations</A>, finite-size
 spheroids and ellipsoids are still treated as 3d particles, rather
 than as disks or ellipses.  This means they have the same moment of
 inertia for a 3d extended object.  When their temperature is
 coomputed, the correct degrees of freedom are used for rotation in a
 2d versus 3d system.
 </P>
 <H5>Pair potentials 
 </H5>
 <P>When a system with extended particles is defined, the particles will
 only rotate and experience torque if the force field computes such
 interactions.  These are the various <A HREF = "pair_style.html">pair styles</A>
 that generate torque:
 </P>
 <UL><LI><A HREF = "pair_gran.html">pair_style gran/history</A>
 <LI><A HREF = "pair_gran.html">pair_style gran/hertzian</A>
 <LI><A HREF = "pair_gran.html">pair_style gran/no_history</A>
 <LI><A HREF = "pair_dipole.html">pair_style dipole/cut</A>
 <LI><A HREF = "pair_gayberne.html">pair_style gayberne</A>
 <LI><A HREF = "pair_resquared.html">pair_style resquared</A>
 <LI><A HREF = "pair_lubricate.html">pair_style lubricate</A> 
 </UL>
 <P>The <A HREF = "pair_gran.html">granular pair styles</A> are used with <A HREF = "atom_style.html">atom_style
 granular</A>.  The <A HREF = "pair_dipole.html">dipole pair style</A>
 is used with <A HREF = "atom_style.html">atom_style dipole</A>.  The
 <A HREF = "pair_gayberne.html">GayBerne</A> and <A HREF = "pair_resquared.html">REsquared</A>
 potentials require particles have a <A HREF = "shape.html">shape</A> and are
 designed for <A HREF = "atom_style.html">ellipsoidal particles</A>.  The
 <A HREF = "pair_lubricate.html">lubrication potential</A> requires that particles
 have a <A HREF = "shape.html">shape</A>.  It can currently only be used with
 extended spherical particles.
 </P>
 <H5>Time integration 
 </H5>
 <P>There are 3 fixes that perform time integration on extended spherical
 particles, meaning the integrators update the rotational orientation
 and angular velocity or angular momentum of the particles:
 </P>
 <UL><LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
 <LI><A HREF = "fix_nvt_sphere.html">fix nvt/sphere</A>
 <LI><A HREF = "fix_npt_sphere.html">fix npt/sphere</A> 
 </UL>
 <P>Likewise, there are 3 fixes that perform time integration on extended
 aspherical particles:
 </P>
 <UL><LI><A HREF = "fix_nve_asphere.html">fix nve/asphere</A>
 <LI><A HREF = "fix_nvt_asphere.html">fix nvt/asphere</A>
 <LI><A HREF = "fix_npt_asphere.html">fix npt/asphere</A> 
 </UL>
 <P>The advantage of these fixes is that those which thermostat the
 particles include the rotational degrees of freedom in the temperature
 calculation and thermostatting.  Other thermostats can be used with
 fix nve/sphere or fix nve/asphere, such as fix langevin or fix
 temp/berendsen, but those thermostats only operate on the
 translational kinetic energy of the extended particles.
 </P>
 <P>Note that for mixtures of point and extended particles, you should
 only use these integration fixes on <A HREF = "group.html">groups</A> which contain
 extended particles.
 </P>
 <H5>Computes, thermodynamics, and dump output 
 </H5>
 <P>There are 4 computes that calculate the temperature or rotational energy
 of extended spherical or aspherical particles:
 </P>
 <UL><LI><A HREF = "compute_temp_sphere.html">compute temp/sphere</A>
 <LI><A HREF = "compute_temp_asphere.html">compute temp/asphere</A>
 <LI><A HREF = "compute_erotate_sphere.html">compute erotate/sphere</A>
 <LI><A HREF = "compute_erotate_asphere.html">compute erotate/asphere</A> 
 </UL>
 <P>These include rotational degrees of freedom in their computation.  If
 you wish the thermodynamic output of temperature or pressure to use
 one of these computes (e.g. for a system entirely composed of extended
 particles), then the compute can be defined and the
 <A HREF = "thermo_modify.html">thermo_modify</A> command used.  Note that by
 default thermodynamic quantities will be calculated with a temperature
 that only includes translational degrees of freedom.  See the
 <A HREF = "thermo_style.html">thermo_style</A> command for details.
 </P>
 <P>The <A HREF = "dump.html">dump custom</A> command can output various attributes of
 extended particles, including the dipole moment (mu), the angular
 velocity (omega), the angular momentum (angmom), the quaternion
 (quat), and the torque (tq) on the particle.
 </P>
 <H5>Rigid bodies composed of extended particles 
 </H5>
 <P>The <A HREF = "fix_rigid.html">fix rigid</A> command treats a collection of
 particles as a rigid body, computes its inertia tensor, sums the total
 force and torque on the rigid body each timestep due to forces on its
 constituent particles, and integrates the motion of the rigid body.
 </P>
 <P>(NOTE: the feature described in the following paragraph has not yet
 been released.  It will be soon.)
 </P>
 <P>If any of the constituent particles of a rigid body are extended
 particles (spheroids or ellipsoids), then their contribution to the
 inertia tensor of the body is different than if they were point
 particles.  This means the rotational dynamics of the rigid body will
 be different.  Thus a model of a dimer is different if the dimer
 consists of two point masses versus two extended sphereoids, even if
 the two particles have the same mass.  Extended particles that
 experience torque due to their interaction with other particles will
 also impart that torque to a rigid body they are part of.
 </P>
 <P>See the "fix rigid" command for example of complex rigid-body models
 it is possible to define in LAMMPS.
 </P>
 <P>Note that the <A HREF = "fix_shake.html">fix shake</A> command can also be used to
 treat 2, 3, or 4 particles as a rigid body, but it always assumes the
 particles are point masses.
 </P>
 <HR>
 
 <A NAME = "4_15"></A><H4>4.15 Output from LAMMPS (thermo, dumps, computes, fixes, variables) 
 </H4>
 <P>There are four basic kinds of LAMMPS output:
 </P>
 <UL><LI><A HREF = "thermo_style.html">Thermodynamic output</A>, which is a list
 of quantities printed every few timesteps to the screen and logfile. 
 
 <LI><A HREF = "dump.html">Dump files</A>, which contain snapshots of atoms and various
 per-atom values and are written at a specified frequency. 
 
 <LI>Certain fixes can output user-specified quantities to files: <A HREF = "fix_ave_time.html">fix
 ave/time</A> for time averaging, <A HREF = "fix_ave_spatial.html">fix
 ave/spatial</A> for spatial averaging, and <A HREF = "fix_print.html">fix
 print</A> for single-line output of
 <A HREF = "variable.html">variables</A>.  Fix print can also output to the
 screen. 
 
 <LI><A HREF = "restart.html">Restart files</A>. 
 </UL>
 <P>A simulation prints one set of thermodynamic output and (optionally)
 restart files.  It can generate any number of dump files and fix
 output files, depending on what <A HREF = "dump.html">dump</A> and <A HREF = "fix.html">fix</A>
 commands you specify.
 </P>
 <P>As discussed below, LAMMPS gives you a variety of ways to determine
 what quantities are computed and printed when the thermodynamics,
 dump, or fix commands listed above perform output.  Throughout this
 discussion, note that users can also <A HREF = "Section_modify.html">add their own computes and fixes
 to LAMMPS</A> which can then generate values that can
 then be output with these commands.
 </P>
 <P>The following sub-sections discuss different LAMMPS command related
 to output and the kind of data they operate on and produce:
 </P>
 <UL><LI><A HREF = "#global">Global/per-atom/local data</A>
 <LI><A HREF = "#scalar">Scalar/vector/array data</A>
 <LI><A HREF = "#thermo">Thermodynamic output</A>
 <LI><A HREF = "#dump">Dump file output</A>
 <LI><A HREF = "#fixoutput">Fixes that write output files</A>
 <LI><A HREF = "#computeoutput">Computes that process output quantities</A>
 <LI><A HREF = "#fixoutput">Fixes that process output quantities</A>
 <LI><A HREF = "#compute">Computes that generate values to output</A>
 <LI><A HREF = "#fix">Fixes that generate values to output</A>
 <LI><A HREF = "#variable">Variables that generate values to output</A>
 <LI><A HREF = "#table">Summary table of output options and data flow between commands</A> 
 </UL>
 <H5><A NAME = "global"></A>Global/per-atom/local data 
 </H5>
 <P>Various output-related commands work with three different styles of
 data: global, per-atom, or local.  A global datum is one or more
 system-wide values, e.g. the temperature of the system.  A per-atom
 datum is one or more values per atom, e.g. the kinetic energy of each
 atom.  Local datums are calculated by each processor based on the
 atoms it owns, but there may be zero or more per atom, e.g. a list of
 bond distances.
 </P>
 <H5><A NAME = "scalar"></A>Scalar/vector/array data 
 </H5>
 <P>Global, per-atom, and local datums can each come in three kinds: a
 single scalar value, a vector of values, or a 2d array of values.  The
 doc page for a "compute" or "fix" or "variable" that generates data
 will specify both the style and kind of data it produces, e.g. a
 per-atom vector.
 </P>
 <P>When a quantity is accessed, as in many of the output commands
 discussed below, it can be referenced via the following bracket
 notation, where ID in this case is the ID of a compute.  The leading
 "c_" would be replaced by "f_" for a fix, or "v_" for a variable:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR><TD >c_ID </TD><TD > entire scalar, vector, or array</TD></TR>
 <TR><TD >c_ID[I] </TD><TD > one element of vector, one column of array</TD></TR>
 <TR><TD >c_ID[I][J] </TD><TD > one element of array 
 </TD></TR></TABLE></DIV>
 
 <P>In other words, using one bracket reduces the dimension of the data
 once (vector -> scalar, array -> vector).  Using two brackets reduces
 the dimension twice (array -> scalar).  Thus a command that uses
 scalar values as input can typically also process elements of a vector
 or array.
 </P>
 <H5><A NAME = "thermo"></A>Thermodynamic output 
 </H5>
 <P>The frequency and format of thermodynamic output is set by the
 <A HREF = "thermo.html">thermo</A>, <A HREF = "thermo_style.html">thermo_style</A>, and
 <A HREF = "thermo_modify.html">thermo_modify</A> commands.  The
 <A HREF = "thermo_style.html">thermo_style</A> command also specifies what values
 are calculated and written out.  Pre-defined keywords can be specified
 (e.g. press, etotal, etc).  Three additional kinds of keywords can
 also be specified (c_ID, f_ID, v_name), where a <A HREF = "compute.html">compute</A>
 or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A> provides the value to be
 output.  In each case, the compute, fix, or variable must generate
 global values for input to the <A HREF = "dump.html">thermo_style custom</A>
 command.
 </P>
 <H5><A NAME = "dump"></A>Dump file output 
 </H5>
 <P>Dump file output is specified by the <A HREF = "dump.html">dump</A> and
 <A HREF = "dump_modify.html">dump_modify</A> commands.  There are several
 pre-defined formats (dump atom, dump xtc, etc).
 </P>
 <P>There is also a <A HREF = "dump.html">dump custom</A> format where the user
 specifies what values are output with each atom.  Pre-defined atom
 attributes can be specified (id, x, fx, etc).  Three additional kinds
 of keywords can also be specified (c_ID, f_ID, v_name), where a
 <A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A>
 provides the values to be output.  In each case, the compute, fix, or
 variable must generate per-atom values for input to the <A HREF = "dump.html">dump
 custom</A> command.
 </P>
 <P>There is also a <A HREF = "dump.html">dump local</A> format where the user specifies
 what local values to output.  A pre-defined index keyword can be
 specified to enumuerate the local values.  Two additional kinds of
 keywords can also be specified (c_ID, f_ID), where a
 <A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A>
 provides the values to be output.  In each case, the compute or fix
 must generate local values for input to the <A HREF = "dump.html">dump local</A>
 command.
 </P>
 <H5><A NAME = "fixoutput"></A>Fixes that write output files 
 </H5>
 <P>Sevarl fixes take various quantities as input and can write output
 files: <A HREF = "fix_ave_time.html">fix ave/time</A>, <A HREF = "fix_ave_spatial.html">fix
 ave/spatial</A>, <A HREF = "fix_ave_histo.html">fix ave/histo</A>,
 <A HREF = "fix_ave_correlate.html">fix ave/correlate</A>, and <A HREF = "fix_print.html">fix
 print</A>.
 </P>
 <P>The <A HREF = "fix_ave_time.html">fix ave/time</A> command enables direct output to
 a file and/or time-averaging of global scalars or vectors.  The user
 specifies one or more quantities as input.  These can be global
 <A HREF = "compute.html">compute</A> values, global <A HREF = "fix.html">fix</A> values, or
 <A HREF = "variable.html">variables</A> of any style except the atom style which
 produces per-atom values.  Since a variable can refer to keywords used
 by the <A HREF = "thermo_style.html">thermo_style custom</A> command (like temp or
 press) and individual per-atom values, a wide variety of quantities
 can be time averaged and/or output in this way.  If the inputs are one
 or more scalar values, then the fix generate a global scalar or vector
 of output.  If the inputs are one or more vector values, then the fix
 generates a global vector or array of output.  The time-averaged
 output of this fix can also be used as input to other output commands.
 </P>
 <P>The <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command enables direct
 output to a file of spatial-averaged per-atom quantities like those
 output in dump files, within 1d layers of the simulation box.  The
 per-atom quantities can be atom density (mass or number) or atom
 attributes such as position, velocity, force.  They can also be
 per-atom quantities calculated by a <A HREF = "compute.html">compute</A>, by a
 <A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>.  The
 spatial-averaged output of this fix can also be used as input to other
 output commands.
 </P>
 <P>The <A HREF = "fix_ave_histo.html">fix ave/histo</A> command enables direct output
 to a file of histogrammed quantities, which can be global or per-atom
 or local quantities.  The histogram output of this fix can also be
 used as input to other output commands.
 </P>
 <P>The <A HREF = "fix_ave_histo.html">fix ave/correlate</A> command enables direct
 output to a file of time-correlated quantities, which can be global
 scalars.  The correlation matrix output of this fix can also be used
 as input to other output commands.
 </P>
 <P>The <A HREF = "fix_print.html">fix print</A> command can generate a line of output
 written to the screen and log file or to a separate file, periodically
 during a running simulation.  The line can contain one or more
 <A HREF = "variable.html">variable</A> values for any style variable except the atom
 style).  As explained above, variables themselves can contain
 references to global values generated by <A HREF = "thermo_style.html">thermodynamic
 keywords</A>, <A HREF = "compute.html">computes</A>,
 <A HREF = "fix.html">fixes</A>, or other <A HREF = "variable.html">variables</A>, or to per-atom
 values for a specific atom.  Thus the <A HREF = "fix_print.html">fix print</A>
 command is a means to output a wide variety of quantities separate
 from normal thermodynamic or dump file output.
 </P>
 <H5><A NAME = "computeoutput"></A>Computes that process output quantities 
 </H5>
 <P>The <A HREF = "compute_reduce.html">compute reduce</A> and <A HREF = "compute_reduce.html">compute
 reduce/region</A> commands take one or more per-atom
 or local vector quantities as inputs and "reduce" them (sum, min, max,
 ave) to scalar quantities.  These are produced as output values which
 can be used as input to other output commands.
 </P>
 <P>The <A HREF = "compute_property_atom.html">compute property/atom</A> command takes a
 list of one or more pre-defined atom attributes (id, x, fx, etc) and
 stores the values in a per-atom vector or array.  These are produced
 as output values which can be used as input to other output commands.
 The list of atom attributes is the same as for the <A HREF = "dump.html">dump
 custom</A> command.
 </P>
 <P>The <A HREF = "compute_property_local.html">compute property/local</A> command takes
 a list of one or more pre-defined local attributes (bond info, angle
 info, etc) and stores the values in a local vector or array.  These
 are produced as output values which can be used as input to other
 output commands.
 </P>
 <P>The <A HREF = "compute_atom_molecule.html">compute atom/molecule</A> command takes a
 list of one or more per-atom quantities (from a compute, fix, per-atom
 variable) and sums the quantities on a per-molecule basis.  It
 produces a global vector or array as output values which can be used
 as input to other output commands.
 </P>
 <H5><A NAME = "fixoutput"></A>Fixes that process output quantities 
 </H5>
 <P>The <A HREF = "fix_ave_atom.html">fix ave/atom</A> command performs time-averaging
 of per-atom vectors.  The per-atom quantities can be atom attributes
 such as position, velocity, force.  They can also be per-atom
 quantities calculated by a <A HREF = "compute.html">compute</A>, by a
 <A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>.  The
 time-averaged per-atom output of this fix can be used as input to
 other output commands.
 </P>
 <P>The <A HREF = "fix_store_state.html">fix store/state</A> command can archive one or
 more per-atom attributes at a particular time, so that the old values
 can be used in a future calculation or output.  The list of atom
 attributes is the same as for the <A HREF = "dump.html">dump custom</A> command,
 including per-atom quantities calculated by a <A HREF = "compute.html">compute</A>,
 by a <A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>.
 The output of this fix can be used as input to other output commands.
 </P>
 <H5><A NAME = "compute"></A>Computes that generate values to output 
 </H5>
 <P>Every <A HREF = "compute.html">compute</A> in LAMMPS produces either global or
 per-atom or local values.  The values can be scalars or vectors or
 arrays of data.  These values can be output using the other commands
 described in this section.  The doc page for each compute command
 describes what it produces.  Computes that produce per-atom or local
 values have the word "atom" or "local" in their style name.  Computes
 without the word "atom" or "local" produce global values.
 </P>
 <H5><A NAME = "fix"></A>Fixes that generate values to output 
 </H5>
 <P>Some <A HREF = "fix.html">fixes</A> in LAMMPS produces either global or per-atom or
 local values which can be accessed by other commands.  The values can
 be scalars or vectors or arrays of data.  These values can be output
 using the other commands described in this section.  The doc page for
 each fix command tells whether it produces any output quantities and
 describes them.
 </P>
 <H5><A NAME = "variable"></A>Variables that generate values to output 
 </H5>
 <P>Every <A HREF = "variable.html">variables</A> defined in an input script generates
 either a global scalar value or a per-atom vector (only atom-style
 variables) when it is accessed.  The formulas used to define equal-
 and atom-style variables can contain references to the thermodynamic
 keywords and to global and per-atom data generated by computes, fixes,
 and other variables.  The values generated by variables can be output
 using the other commands described in this section.
 </P>
 <H5><A NAME = "table"></A>Summary table of output options and data flow between commands 
 </H5>
 <P>This table summarizes the various commands that can be used for
 generating output from LAMMPS.  Each command produces output data of
 some kind and/or writes data to a file.  Most of the commands can take
 data from other commands as input.  Thus you can link many of these
 commands together in pipeline form, where data produced by one command
 is used as input to another command and eventually written to the
 screen or to a file.  Note that to hook two commands together the
 output and input data types must match, e.g. global/per-atom/local
 data and scalar/vector/array data.
 </P>
 <P>Also note that, as described above, when a command takes a scalar as
 input, that could be an element of a vector or array.  Likewise a
 vector input could be a column of an array.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR><TD >Command</TD><TD > Input</TD><TD > Output</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "thermo_style.html">thermo_style custom</A></TD><TD > global scalars</TD><TD > screen, log file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "dump.html">dump custom</A></TD><TD > per-atom vectors</TD><TD > dump file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "dump.html">dump local</A></TD><TD > local vectors</TD><TD > dump file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_print.html">fix print</A></TD><TD > global scalar from variable</TD><TD > screen, file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "print.html">print</A></TD><TD > global scalar from variable</TD><TD > screen</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "compute.html">computes</A></TD><TD > N/A</TD><TD > global/per-atom/local scalar/vector/array</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix.html">fixes</A></TD><TD > N/A</TD><TD > global/per-atom/local scalar/vector/array</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "variable.html">variables</A></TD><TD > global scalars, per-atom vectors</TD><TD > global scalar, per-atom vector</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "compute_reduce.html">compute reduce</A></TD><TD > per-atom/local vectors</TD><TD > global scalar/vector</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "compute_property_atom.html">compute property/atom</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "compute_property_local.html">compute property/local</A></TD><TD > local vectors</TD><TD > local vector/array</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "compute_atom_molecule.html">compute atom/molecule</A></TD><TD > per-atom vectors</TD><TD > global vector/array</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_ave_atom.html">fix ave/atom</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_ave_time.html">fix ave/time</A></TD><TD > global scalars/vectors</TD><TD > global scalar/vector/array, file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_ave_spatial.html">fix ave/spatial</A></TD><TD > per-atom vectors</TD><TD > global array, file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_ave_histo.html">fix ave/histo</A></TD><TD > global/per-atom/local scalars and vectors</TD><TD > global array, file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_ave_correlate.html">fix ave/correlate</A></TD><TD > global scalars</TD><TD > global array, file</TD><TD ></TD></TR>
 <TR><TD ><A HREF = "fix_store_state.html">fix store/state</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
 <TR><TD >
 </TD></TR></TABLE></DIV>
 
 <HR>
 
 <A NAME = "4_16"></A><H4>4.16 Thermostatting, barostatting, and computing temperature 
 </H4>
 <P>Thermostatting means controlling the temperature of particles in an MD
 simulation.  Barostatting means controlling the pressure.  Since the
 pressure includes a kinetic component due to particle velocities, both
 these operations require calculation of the temperature.  Typically a
 target temperature (T) and/or pressure (P) is specified by the user,
 and the thermostat or barostat attempts to equilibrate the system to
 the requested T and/or P.
 </P>
 <P>Temperature is computed as kinetic energy divided by some number of
 degrees of freedom (and the Boltzmann constant).  Since kinetic energy
 is a function of particle velocity, there is often a need to
 distinguish between a particle's advection velocity (due to some
 aggregate motiion of particles) and its thermal velocity.  The sum of
 the two is the particle's total velocity, but the latter is often what
 is wanted to compute a temperature.
 </P>
 <P>LAMMPS has several options for computing temperatures, any of which
 can be used in thermostatting and barostatting.  These <A HREF = "compute.html">compute
 commands</A> calculate temperature, and the <A HREF = "compute_pressure.html">compute
 pressure</A> command calculates pressure.
 </P>
 <UL><LI><A HREF = "compute_temp.html">compute temp</A>
 <LI><A HREF = "compute_temp_sphere.html">compute temp/sphere</A>
 <LI><A HREF = "compute_temp_asphere.html">compute temp/asphere</A>
 <LI><A HREF = "compute_temp_com.html">compute temp/com</A>
 <LI><A HREF = "compute_temp_deform.html">compute temp/deform</A>
 <LI><A HREF = "compute_temp_partial.html">compute temp/partial</A>
 <LI><A HREF = "compute_temp_profile.html">compute temp/profile</A>
 <LI><A HREF = "compute_temp_ramp.html">compute temp/ramp</A>
 <LI><A HREF = "compute_temp_region.html">compute temp/region</A> 
 </UL>
 <P>All but the first 3 calculate velocity biases (i.e. advection
 velocities) that are removed when computing the thermal temperature.
 <A HREF = "compute_temp_sphere.html">Compute temp/sphere</A> and <A HREF = "compute_temp_asphere.html">compute
 temp/asphere</A> compute kinetic energy for
 extended particles that includes rotational degrees of freedom.  They
 both allow, as an extra argument, which is another temperature compute
 that subtracts a velocity bias.  This allows the translational
 velocity of extended spherical or aspherical particles to be adjusted
 in prescribed ways.
 </P>
 <P>Thermostatting in LAMMPS is performed by <A HREF = "fix.html">fixes</A>, or in one
 case by a pair style.  Four thermostatting fixes are currently
 available: Nose-Hoover (nvt), Berendsen, Langevin, and direct
 rescaling (temp/rescale).  Dissipative particle dynamics (DPD)
 thermostatting can be invoked via the <I>dpd/tstat</I> pair style:
 </P>
 <UL><LI><A HREF = "fix_nh.html">fix nvt</A>
 <LI><A HREF = "fix_nvt_sphere.html">fix nvt/sphere</A>
 <LI><A HREF = "fix_nvt_asphere.html">fix nvt/asphere</A>
 <LI><A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A>
 <LI><A HREF = "fix_temp_berendsen.html">fix temp/berendsen</A>
 <LI><A HREF = "fix_langevin.html">fix langevin</A>
 <LI><A HREF = "fix_temp_rescale.html">fix temp/rescale</A>
 <LI><A HREF = "pair_dpd.html">pair_style dpd/tstat</A> 
 </UL>
 <P><A HREF = "fix_nh.html">Fix nvt</A> only thermostats the translational velocity of
 particles.  <A HREF = "fix_nvt_sllod.html">Fix nvt/sllod</A> also does this, except
 that it subtracts out a velocity bias due to a deforming box and
 integrates the SLLOD equations of motion.  See the <A HREF = "#4_13">NEMD
 simulations</A> section of this page for further details.  <A HREF = "fix_nvt_sphere.html">Fix
 nvt/sphere</A> and <A HREF = "fix_nvt_asphere.html">fix
 nvt/asphere</A> thermostat not only translation
 velocities but also rotational velocities for spherical and aspherical
 particles.
 </P>
 <P>DPD thermostatting alters pairwise interactions in a manner analagous
 to the per-particle thermostatting of <A HREF = "fix_langevin.html">fix
 langevin</A>.
 </P>
 <P>Any of the thermostatting fixes can use temperature computes that
 remove bias for two purposes: (a) computing the current temperature to
 compare to the requested target temperature, and (b) adjusting only
 the thermal temperature component of the particle's velocities.  See
 the doc pages for the individual fixes and for the
 <A HREF = "fix_modify.html">fix_modify</A> command for instructions on how to assign
 a temperature compute to a thermostatting fix.  For example, you can
 apply a thermostat to only the x and z components of velocity by using
 it in conjunction with <A HREF = "compute_temp_partial.html">compute
 temp/partial</A>.
 </P>
 <P>IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning
 they update the velocities and positions of particles due to forces
 and velocities respectively.  The other thermostat fixes only adjust
 velocities; they do NOT perform time integration updates.  Thus they
 should be used in conjunction with a constant NVE integration fix such
 as these:
 </P>
 <UL><LI><A HREF = "fix_nve.html">fix nve</A>
 <LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
 <LI><A HREF = "fix_nve_asphere.html">fix nve/asphere</A> 
 </UL>
 <P>Barostatting in LAMMPS is also performed by <A HREF = "fix.html">fixes</A>.  Two
 barosttating methods are currently available: Nose-Hoover (npt and
 nph) and Berendsen:
 </P>
 <UL><LI><A HREF = "fix_nh.html">fix npt</A>
 <LI><A HREF = "fix_npt_sphere.html">fix npt/sphere</A>
 <LI><A HREF = "fix_npt_asphere.html">fix npt/asphere</A>
 <LI><A HREF = "fix_nh.html">fix nph</A>
 <LI><A HREF = "fix_press_berendsen.html">fix press/berendsen</A> 
 </UL>
 <P>The <A HREF = "fix_nh.html">fix npt</A> commands include a Nose-Hoover thermostat
 and barostat.  <A HREF = "fix_nh.html">Fix nph</A> is just a Nose/Hoover barostat;
 it does no thermostatting.  Both <A HREF = "fix_nh.html">fix nph</A> and <A HREF = "fix_press_berendsen.html">fix
 press/bernendsen</A> can be used in conjunction
 with any of the thermostatting fixes.
 </P>
 <P>As with the thermostats, <A HREF = "fix_nh.html">fix npt</A> and <A HREF = "fix_nh.html">fix
 nph</A> only use translational motion of the particles in
 computing T and P and performing thermo/barostatting.  <A HREF = "fix_npt_sphere.html">Fix
 npt/sphere</A> and <A HREF = "fix_npt_asphere.html">fix
 npt/asphere</A> thermo/barostat using not only
 translation velocities but also rotational velocities for spherical
 and aspherical particles.
 </P>
 <P>All of the barostatting fixes use the <A HREF = "compute_pressure.html">compute
 pressure</A> compute to calculate a current
 pressure.  By default, this compute is created with a simple <A HREF = "compute_temp.html">compute
 temp</A> (see the last argument of the <A HREF = "compute_pressure.html">compute
 pressure</A> command), which is used to calculated
 the kinetic componenet of the pressure.  The barostatting fixes can
 also use temperature computes that remove bias for the purpose of
 computing the kinetic componenet which contributes to the current
 pressure.  See the doc pages for the individual fixes and for the
 <A HREF = "fix_modify.html">fix_modify</A> command for instructions on how to assign
 a temperature or pressure compute to a barostatting fix.
 </P>
 <P>IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods (<A HREF = "fix_nh.html">fix
 npt</A> and <A HREF = "fix_nh.html">fix nph</A>) perform time
 integration.  <A HREF = "fix_press_berendsen.html">Fix press/berendsen</A> does NOT,
 so it should be used with one of the constant NVE fixes or with one of
 the NVT fixes.
 </P>
 <P>Finally, thermodynamic output, which can be setup via the
 <A HREF = "thermo_style.html">thermo_style</A> command, often includes temperature
 and pressure values.  As explained on the doc page for the
 <A HREF = "thermo_style.html">thermo_style</A> command, the default T and P are
 setup by the thermo command itself.  They are NOT the ones associated
 with any thermostatting or barostatting fix you have defined or with
 any compute that calculates a temperature or pressure.  Thus if you
 want to view these values of T and P, you need to specify them
 explicitly via a <A HREF = "thermo_style.html">thermo_style custom</A> command.  Or
 you can use the <A HREF = "thermo_modify.html">thermo_modify</A> command to
 re-define what temperature or pressure compute is used for default
 thermodynamic output.
 </P>
 <HR>
 
 <A NAME = "4_17"></A><H4>4.17 Walls 
 </H4>
 <P>Walls in an MD simulation are typically used to bound particle motion,
 i.e. to serve as a boundary condition.
 </P>
 <P>Walls in LAMMPS can be of rough (made of particles) or idealized
 surfaces.  Ideal walls can be smooth, generating forces only in the
 normal direction, or frictional, generating forces also in the
 tangential direction.
 </P>
 <P>Rough walls, built of particles, can be created in various ways.  The
 particles themselves can be generated like any other particle, via the
 <A HREF = "lattice.html">lattice</A> and <A HREF = "create_atoms.html">create_atoms</A> commands,
 or read in via the <A HREF = "read_data.html">read_data</A> command.
 </P>
 <P>Their motion can be constrained by many different commands, so that
 they do not move at all, move together as a group at constant velocity
 or in response to a net force acting on them, move in a prescribed
 fashion (e.g. rotate around a point), etc.  Note that if a time
 integration fix like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix nvt</A>
 is not used with the group that contains wall particles, their
 positions and velocities will not be updated.
 </P>
 <UL><LI><A HREF = "fix_aveforce.html">fix aveforce</A> - set force on particles to average value, so they move together
 <LI><A HREF = "fix_setforce.html">fix setforce</A> - set force on particles to a value, e.g. 0.0
 <LI><A HREF = "fix_freeze.html">fix freeze</A> - freeze particles for use as granular walls
 <LI><A HREF = "fix_nve_noforce.html">fix nve/noforce</A> - advect particles by their velocity, but without force
 <LI><A HREF = "fix_move.html">fix move</A> - prescribe motion of particles by a linear velocity, oscillation, rotation, variable 
 </UL>
 <P>The <A HREF = "fix_move.html">fix move</A> command offers the most generality, since
 the motion of individual particles can be specified with
 <A HREF = "variable.html">variable</A> formula which depends on time and/or the
 particle position.
 </P>
 <P>For rough walls, it may be useful to turn off pairwise interactions
 between wall particles via the <A HREF = "neigh_modify.html">neigh_modify
 exclude</A> command.
 </P>
 <P>Rough walls can also be created by specifying frozen particles that do
 not move and do not interact with mobile particles, and then tethering
 other particles to the fixed particles, via a <A HREF = "bond_style.html">bond</A>.
 The bonded particles do interact with other mobile particles.
 </P>
 <P>Idealized walls can be specified via several fix commands.  <A HREF = "fix_wall_gran.html">Fix
 wall/gran</A> creates frictional walls for use with
 granular particles; all the other commands create smooth walls.
 </P>
 <UL><LI><A HREF = "fix_wall_reflect.html">fix wall/reflect</A> - reflective flat walls
 <LI><A HREF = "fix_wall.html">fix wall/lj93</A> - flat walls, with Lennard-Jones 9/3 potential
 <LI><A HREF = "fix_wall.html">fix wall/lj126</A> - flat walls, with Lennard-Jones 12/6 potential
 <LI><A HREF = "fix_wall.html">fix wall/colloid</A> - flat walls, with <A HREF = "pair_colloid.html">pair_style colloid</A> potential
 <LI><A HREF = "fix_wall.html">fix wall/harmonic</A> - flat walls, with repulsive harmonic spring potential
 <LI><A HREF = "fix_wall_region.html">fix wall/region</A> - use region surface as wall
 <LI><A HREF = "fix_wall_gran.html">fix wall/gran</A> - flat or curved walls with <A HREF = "pair_gran.html">pair_style granular</A> potential 
 </UL>
 <P>The <I>lj93</I>, <I>lj126</I>, <I>colloid</I>, and <I>harmonic</I> styles all allow the
 flat walls to move with a constant velocity, or oscillate in time.
 The <A HREF = "fix_wall_region.html">fix wall/region</A> command offers the most
 generality, since the region surface is treated as a wall, and the
 geometry of the region can be a simple primitive volume (e.g. a
 sphere, or cube, or plane), or a complex volume made from the union
 and intersection of primitive volumes.  <A HREF = "region.html">Regions</A> can also
 specify a volume "interior" or "exterior" to the specified primitive
 shape or <I>union</I> or <I>intersection</I>.  <A HREF = "region.html">Regions</A> can also be
 "dynamic" meaning they move with constant velocity, oscillate, or
 rotate.
 </P>
 <P>The only frictional idealized walls currently in LAMMPS are flat or
 curved surfaces specified by the <A HREF = "fix_wall_gran.html">fix wall/gran</A>
 command.  At some point we plan to allow regoin surfaces to be used as
 frictional walls, as well as triangulated surfaces.
 </P>
 <HR>
 
 <A NAME = "4_18"></A><H4>4.18 Elastic constants 
 </H4>
 <P>Elastic constants characterize the stiffness of a material. The formal
 definition is provided by the linear relation that holds between the
 stress and strain tensors in the limit of infinitesimal deformation.
 In tensor notation, this is expressed as s_ij = C_ijkl * e_kl, where
 the repeated indices imply summation. s_ij are the elements of the
 symmetric stress tensor. e_kl are the elements of the symmetric strain
 tensor. C_ijkl are the elements of the fourth rank tensor of elastic
 constants. In three dimensions, this tensor has 3^4=81 elements. Using
 Voigt notation, the tensor can be written as a 6x6 matrix, where C_ij
 is now the derivative of s_i w.r.t. e_j. Because s_i is itself a
 derivative w.r.t. e_i, it follows that C_ij is also symmetric, with at
 most 7*6/2 = 21 distinct elements.
 </P>
 <P>At zero temperature, it is easy to estimate these derivatives by
 deforming the cell in one of the six directions using the command
 <A HREF = "displace_box.html">displace_box</A> and measuring the change in the
 stress tensor. A general-purpose script that does this is given in the
 examples/elastic directory described in <A HREF = "Section_example.html">this
 section</A>.
 </P>
 <P>Calculating elastic constants at finite temperature is more
 challenging, because it is necessary to run a simulation that perfoms
 time averages of differential properties. One way to do this is to
 measure the change in average stress tensor in an NVT simulations when
 the cell volume undergoes a finite deformation. In order to balance
 the systematic and statistical errors in this method, the magnitude of
 the deformation must be chosen judiciously, and care must be taken to
 fully equilibrate the deformed cell before sampling the stress
 tensor. Another approach is to sample the triclinic cell fluctuations
 that occur in an NPT simulation. This method can also be slow to
 converge and requires careful post-processing <A HREF = "#Shinoda">(Shinoda)</A>
 </P>
 <HR>
 
-<A NAME = "4_19"></A><H4>4.19 Computing free energies from thermodynamic integration 
-</H4>
-<P>Thermodynamic integration is a widely used method to compute free
-energies from atomistic simulations.  LAMMPS can be used to run
-thermodynamic integration calculations using the methods discussed in
-this section and the <A HREF = "fix_adapt.html">fix adapt</A> command.  Currently,
-it is capable of the transformations essential for computing melting
-points using the pseudo-supercritical path method developed by <A HREF = "#Eike_Maginn">Eike
-and Maginn</A>.
-</P>
-<P>See the examples/TI directory for more information and sample files
-that compute a melting point using the techniques described in this
-section.  That directory has its own README file.  See also the paper
-by <A HREF = "#Jayaraman">Jayaraman</A> for an example of using this implementation
-of thermodynamic integration in LAMMPS to compute melting points of
-alkali nitrate salts, using the steps outlined here.
-</P>
-<P>In this method, three intermediate "pseudo-supercritical" states are
-accessed in the transformation between the liquid and solid
-phases. These pseudo-states are a weak liquid, a dense weak liquid,
-and an ordered weak phase. The transformation between the liquid and
-solid states can also be driven uisng the <A HREF = "fix_adapt.html">fix adapt</A>
-command.
-</P>
-<P>For the transformation from the liquid to the weak phase, the
-intermolecular interactions need to be weakened.  Appropriate scale
-factors, computed by variables you define, and applied to pair styles
-by <A HREF = "fix_adapt.html">fix adapt</A>, can be used to do this, as in the
-example scripts.  The <A HREF = "compute_ti.html">compute ti</A> command can
-accumulate the value of dU/d<I>lambda</I>.  See <A HREF = "#Jayaraman_Maginn">Jayaraman and
-Maginn</A> for more information about calculating a
-free energy from dU/d<I>lambda</I>.
-</P>
-<P>IMPORTANT NOTE: The pair styles that fix adapt can scale on-the-fly
-are listed on the <A HREF = "fix_adapt">fix adapt</A> doc page.  interaction scaling
-is desired.  If a pair style is not on that list, it is generally
-quite easy to add an extract() method to the pair style, to enable fix
-adapt to rescale it.
-</P>
-<P>Step 2 is the transformation of the simulation box density from the
-liquid phase to that of the equilibrated crystal.  The parameters for
-box1 and box2 should be obtained from equilibrated NPT simulations of
-the liquid and crystal phases and used in a <A HREF = "fix_deform.html">fix
-deform</A> command to change the box size and/or shape.
-It also advisable to use <A HREF = "fix_adapt.html">fix adapt</A> on the pair styles
-to prevent overlaps which may occur during the box transformation.
-</P>
-<P>In step 3, the dense, weak system is transformed to an ordered state,
-which has the same ordering as in the equilibrated crystal.  Ordering
-is achieved by introducing an attractive potential between atoms and
-lattice sites.  These lattice sites can be calculated as the mean
-positions of the atoms in an equilibrium simulation of the
-crystal. The <A HREF = "pair_gauss.html">pair/gauss</A> command can be used to
-introduce an attractive Gaussian potential between the atoms and their
-corresponding lattice sites. The prefactor of the Gaussian pair
-potential can be scaled by <A HREF = "fix_adapt.html">fix adapt</A> to turn on the
-attractions.  Again, the quantity dU/d<I>lambda</I> can be tracked via the
-<A HREF = "compute_ti.html">compute ti</A> command.
-</P>
-<P>Step 4 is the transformation of the ordered state to the final
-crystal.  In this step, the intermolecular interactions are scaled
-back to full strength, while the Gaussian tethers are removed, all via
-<A HREF = "fix_adapt.html">fix adapt</A>.
-</P>
-<HR>
-
 <HR>
 
 <A NAME = "Berendsen"></A>
 
 <P><B>(Berendsen)</B> Berendsen, Grigera, Straatsma, J Phys Chem, 91,
 6269-6271 (1987).
 </P>
 <A NAME = "Cornell"></A>
 
 <P><B>(Cornell)</B> Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
 Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
 </P>
 <A NAME = "Horn"></A>
 
 <P><B>(Horn)</B> Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
 J Chem Phys, 120, 9665 (2004).
 </P>
 <A NAME = "MacKerell"></A>
 
 <P><B>(MacKerell)</B> MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
 Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
 </P>
+<A NAME = "Mayo"></A>
+
+<P><B>(Mayo)</B> Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
+(1990).
+</P>
 <A NAME = "Jorgensen"></A>
 
 <P><B>(Jorgensen)</B> Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
 Phys, 79, 926 (1983).
 </P>
 <A NAME = "Price"></A>
 
 <P><B>(Price)</B> Price and Brooks, J Chem Phys, 121, 10096 (2004).
 </P>
 <A NAME = "Shinoda"></A>
 
 <P><B>(Shinoda)</B> Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
 </P>
-<A NAME = "Eike_Maginn"></A>
-
-<P><B>(Eike and Maginn)</B> Eike and Maginn, J Chem Phys, 124,
-164503 (2006). 
-</P>
-<A NAME = "Jayaraman_Maginn"></A>
-
-<P><B>(Jayaraman and Maginn)</B> Jayaraman and Maginn, Journal of Chemical Physics, 
-127, 214504 (2007).
-</P>
-<A NAME = "Jayaraman"></A>
-
-<P><B>(Jayaraman)</B> Jayaraman, Thompson, von Lilienfeld and Maginn, Industrial 
-and Engineering Chemistry Research, 49, 559-571 (2010).   
-</P>
 </HTML>
diff --git a/doc/Section_howto.txt b/doc/Section_howto.txt
index 22ab363dc..1be717ebd 100644
--- a/doc/Section_howto.txt
+++ b/doc/Section_howto.txt
@@ -1,1704 +1,1668 @@
 "Previous Section"_Section_commands.html - "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc - "Next Section"_Section_example.html :c
 
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 
 :line 
 
 4. How-to discussions :h3
 
 The following sections describe what commands can be used to perform
 certain kinds of LAMMPS simulations.
 
 4.1 "Restarting a simulation"_#4_1
 4.2 "2d simulations"_#4_2
-4.3 "CHARMM and AMBER force fields"_#4_3
+4.3 "CHARMM, AMBER, and DREIDING force fields"_#4_3
 4.4 "Running multiple simulations from one input script"_#4_4
 4.5 "Multi-replica simulations"_#4_5
 4.6 "Granular models"_#4_6
 4.7 "TIP3P water model"_#4_7
 4.8 "TIP4P water model"_#4_8
 4.9 "SPC water model"_#4_9
 4.10 "Coupling LAMMPS to other codes"_#4_10
 4.11 "Visualizing LAMMPS snapshots"_#4_11
 4.12 "Triclinic (non-orthogonal) simulation boxes"_#4_12
 4.13 "NEMD simulations"_#4_13
 4.14 "Extended spherical and aspherical particles"_#4_14
 4.15 "Output from LAMMPS (thermo, dumps, computes, fixes, variables)"_#4_15
 4.16 "Thermostatting, barostatting and computing temperature"_#4_16
 4.17 "Walls"_#4_17
-4.18 "Elastic constants"_#4_18
-4.19 "Computing free energies from thermodyanmic integration"_#4_19 :all(b)
+4.18 "Elastic constants"_#4_18 :all(b)
 
 The example input scripts included in the LAMMPS distribution and
 highlighted in "this section"_Section_example.html also show how to
 setup and run various kinds of problems.
 
 :line
 
 4.1 Restarting a simulation :link(4_1),h4
 
 There are 3 ways to continue a long LAMMPS simulation.  Multiple
 "run"_run.html commands can be used in the same input script.  Each
 run will continue from where the previous run left off.  Or binary
 restart files can be saved to disk using the "restart"_restart.html
 command.  At a later time, these binary files can be read via a
 "read_restart"_read_restart.html command in a new script.  Or they can
 be converted to text data files and read by a
 "read_data"_read_data.html command in a new script.  "This
 section"_Section_tools.html discusses the {restart2data} tool that is
 used to perform the conversion.
 
 Here we give examples of 2 scripts that read either a binary restart
 file or a converted data file and then issue a new run command to
 continue where the previous run left off.  They illustrate what
 settings must be made in the new script.  Details are discussed in the
 documentation for the "read_restart"_read_restart.html and
 "read_data"_read_data.html commands.
 
 Look at the {in.chain} input script provided in the {bench} directory
 of the LAMMPS distribution to see the original script that these 2
 scripts are based on.  If that script had the line
 
 restart	        50 tmp.restart :pre
 
 added to it, it would produce 2 binary restart files (tmp.restart.50
 and tmp.restart.100) as it ran.
 
 This script could be used to read the 1st restart file and re-run the
 last 50 timesteps:
 
 read_restart	tmp.restart.50 :pre
 
 neighbor	0.4 bin
 neigh_modify	every 1 delay 1 :pre
 
 fix		1 all nve
 fix		2 all langevin 1.0 1.0 10.0 904297 :pre
 
 timestep	0.012 :pre
 
 run		50 :pre
 
 Note that the following commands do not need to be repeated because
 their settings are included in the restart file: {units, atom_style,
 special_bonds, pair_style, bond_style}.  However these commands do
 need to be used, since their settings are not in the restart file:
 {neighbor, fix, timestep}.
 
 If you actually use this script to perform a restarted run, you will
 notice that the thermodynamic data match at step 50 (if you also put a
 "thermo 50" command in the original script), but do not match at step
 100.  This is because the "fix langevin"_fix_langevin.html command
 uses random numbers in a way that does not allow for perfect restarts.
 
 As an alternate approach, the restart file could be converted to a data
 file using this tool:
 
 restart2data tmp.restart.50 tmp.restart.data :pre
 
 Then, this script could be used to re-run the last 50 steps:
 
 units		lj
 atom_style	bond
 pair_style	lj/cut 1.12
 pair_modify	shift yes
 bond_style	fene
 special_bonds   0.0 1.0 1.0 :pre
 
 read_data	tmp.restart.data :pre
 
 neighbor	0.4 bin
 neigh_modify	every 1 delay 1 :pre
 
 fix		1 all nve
 fix		2 all langevin 1.0 1.0 10.0 904297 :pre
 
 timestep	0.012 :pre
 
 reset_timestep	50
 run		50 :pre
 
 Note that nearly all the settings specified in the original {in.chain}
 script must be repeated, except the {pair_coeff} and {bond_coeff}
 commands since the new data file lists the force field coefficients.
 Also, the "reset_timestep"_reset_timestep.html command is used to tell
 LAMMPS the current timestep.  This value is stored in restart files,
 but not in data files.
 
 :line
 
 4.2 2d simulations :link(4_2),h4
 
 Use the "dimension"_dimension.html command to specify a 2d simulation.
 
 Make the simulation box periodic in z via the "boundary"_boundary.html
 command.  This is the default.
 
 If using the "create box"_create_box.html command to define a
 simulation box, set the z dimensions narrow, but finite, so that the
 create_atoms command will tile the 3d simulation box with a single z
 plane of atoms - e.g.
 
 "create box"_create_box.html 1 -10 10 -10 10 -0.25 0.25 :pre
 
 If using the "read data"_read_data.html command to read in a file of
 atom coordinates, set the "zlo zhi" values to be finite but narrow,
 similar to the create_box command settings just described.  For each
 atom in the file, assign a z coordinate so it falls inside the
 z-boundaries of the box - e.g. 0.0.
 
 Use the "fix enforce2d"_fix_enforce2d.html command as the last
 defined fix to insure that the z-components of velocities and forces
 are zeroed out every timestep.  The reason to make it the last fix is
 so that any forces induced by other fixes will be zeroed out.
 
 Many of the example input scripts included in the LAMMPS distribution
 are for 2d models.
 
 IMPORTANT NOTE: Some models in LAMMPS treat particles as extended
 spheres, as opposed to point particles.  In 2d, the particles will
 still be spheres, not disks, meaning their moment of inertia will be
 the same as in 3d.
 
 :line
 
-4.3 CHARMM and AMBER force fields :link(4_3),h4
+4.3 CHARMM, AMBER, and DREIDING force fields :link(4_3),h4
 
-There are many different ways to compute forces in the "CHARMM"_charmm
-and "AMBER"_amber molecular dynamics codes, only some of which are
-available as options in LAMMPS.  A force field has 2 parts: the
-formulas that define it and the coefficients used for a particular
-system.  Here we only discuss formulas implemented in LAMMPS.  Setting
+A force field has 2 parts: the formulas that define it and the
+coefficients used for a particular system.  Here we only discuss
+formulas implemented in LAMMPS that correspond to formulas commonly
+used in the CHARMM, AMBER, and DREIDING force fields.  Setting
 coefficients is done in the input data file via the
 "read_data"_read_data.html command or in the input script with
 commands like "pair_coeff"_pair_coeff.html or
-"bond_coeff"_bond_coeff.html.  See "this section"_Section_tools.html for
-additional tools that can use CHARMM or AMBER to assign force field
-coefficients and convert their output into LAMMPS input.
+"bond_coeff"_bond_coeff.html.  See "this section"_Section_tools.html
+for additional tools that can use CHARMM or AMBER to assign force
+field coefficients and convert their output into LAMMPS input.
 
 See "(MacKerell)"_#MacKerell for a description of the CHARMM force
 field.  See "(Cornell)"_#Cornell for a description of the AMBER force
 field.
 
 :link(charmm,http://www.scripps.edu/brooks)
 :link(amber,http://amber.scripps.edu)
 
 These style choices compute force field formulas that are consistent
 with common options in CHARMM or AMBER.  See each command's
 documentation for the formula it computes.
 
-"bond_style"_bond_style.html harmonic
-"angle_style"_angle_style.html charmm
-"dihedral_style"_dihedral_style.html charmm
-"pair_style"_pair_style.html lj/charmm/coul/charmm
-"pair_style"_pair_style.html lj/charmm/coul/charmm/implicit
-"pair_style"_pair_style.html lj/charmm/coul/long :ul
+"bond_style"_bond_harmonic.html harmonic
+"angle_style"_angle_charmm.html charmm
+"dihedral_style"_dihedral_charmm.html charmm
+"pair_style"_pair_charmm.html lj/charmm/coul/charmm
+"pair_style"_pair_charmm.html lj/charmm/coul/charmm/implicit
+"pair_style"_pair_charmm.html lj/charmm/coul/long :ul
 
 "special_bonds"_special_bonds.html charmm
 "special_bonds"_special_bonds.html amber :ul
 
+DREIDING is a generic force field developed by the "Goddard
+group"_http://www.wag.caltech.edu at Caltech and is useful for
+predicting structures and dynamics of organic, biological and
+main-group inorganic molecules. The philosophy in DREIDING is to use
+general force constants and geometry parameters based on simple
+hybridization considerations, rather than individual force constants
+and geometric parameters that depend on the particular combinations of
+atoms involved in the bond, angle, or torsion terms. DREIDING has an
+"explicit hydrogen bond term"_pair_hbond_dreiding.html to describe
+interactions involving a hydrogen atom (H___A) on very electronegative
+atoms (N, O, F).
+
+See "(Mayo)"_#Mayo for a description of the DREIDING force field
+
+These style choices compute force field formulas that are consistent
+with the DREIDING force field.  See each command's
+documentation for the formula it computes.
+
+"bond_style"_bond_harmonic.html harmonic
+"bond_style"_bond_morse.html morse :ul
+
+"angle_style"_angle_harmonic.html harmonic
+"angle_style"_angle_cosine.html cosine
+"angle_style"_angle_cosine_periodic.html cosine/periodic :ul
+
+"dihedral_style"_dihedral_charmm.html charmm
+"improper_style"_improper_umbrella.html umbrella :ul
+
+"pair_style"_pair_buck.html buck
+"pair_style"_pair_buck.html buck/coul/cut
+"pair_style"_pair_buck.html buck/coul/long
+"pair_style"_pair_lj.html lj/cut
+"pair_style"_pair_lj.html lj/cut/coul/cut
+"pair_style"_pair_lj.html lj/cut/coul/long :ul
+
+"pair_style"_pair_hbond_dreiding.html hbond/dreiding/lj
+"pair_style"_pair_hbond_dreiding.html hbond/dreiding/morse :ul
+
+"special_bonds"_special_bonds.html dreiding :ul
+
 :line
 
 4.4 Running multiple simulations from one input script :link(4_4),h4
 
 This can be done in several ways.  See the documentation for
 individual commands for more details on how these examples work.
 
 If "multiple simulations" means continue a previous simulation for
 more timesteps, then you simply use the "run"_run.html command
 multiple times.  For example, this script
 
 units lj
 atom_style atomic
 read_data data.lj
 run 10000
 run 10000
 run 10000
 run 10000
 run 10000 :pre
 
 would run 5 successive simulations of the same system for a total of
 50,000 timesteps.
 
 If you wish to run totally different simulations, one after the other,
 the "clear"_clear.html command can be used in between them to
 re-initialize LAMMPS.  For example, this script
 
 units lj
 atom_style atomic
 read_data data.lj
 run 10000
 clear
 units lj
 atom_style atomic
 read_data data.lj.new
 run 10000 :pre
 
 would run 2 independent simulations, one after the other.
 
 For large numbers of independent simulations, you can use
 "variables"_variable.html and the "next"_next.html and
 "jump"_jump.html commands to loop over the same input script
 multiple times with different settings.  For example, this
 script, named in.polymer
 
 variable d index run1 run2 run3 run4 run5 run6 run7 run8
 shell cd $d
 read_data data.polymer
 run 10000
 shell cd ..
 clear
 next d
 jump in.polymer :pre
 
 would run 8 simulations in different directories, using a data.polymer
 file in each directory.  The same concept could be used to run the
 same system at 8 different temperatures, using a temperature variable
 and storing the output in different log and dump files, for example
 
 variable a loop 8
 variable t index 0.8 0.85 0.9 0.95 1.0 1.05 1.1 1.15
 log log.$a
 read data.polymer
 velocity all create $t 352839
 fix 1 all nvt $t $t 100.0
 dump 1 all atom 1000 dump.$a
 run 100000
 next t
 next a
 jump in.polymer :pre
 
 All of the above examples work whether you are running on 1 or
 multiple processors, but assumed you are running LAMMPS on a single
 partition of processors.  LAMMPS can be run on multiple partitions via
 the "-partition" command-line switch as described in "this
 section"_Section_start.html#2_6 of the manual.
 
 In the last 2 examples, if LAMMPS were run on 3 partitions, the same
 scripts could be used if the "index" and "loop" variables were
 replaced with {universe}-style variables, as described in the
 "variable"_variable.html command.  Also, the "next t" and "next a"
 commands would need to be replaced with a single "next a t" command.
 With these modifications, the 8 simulations of each script would run
 on the 3 partitions one after the other until all were finished.
 Initially, 3 simulations would be started simultaneously, one on each
 partition.  When one finished, that partition would then start
 the 4th simulation, and so forth, until all 8 were completed.
 
 :line
 
 4.5 Multi-replica simulations :link(4_5),h4
 
 Several commands in LAMMPS run mutli-replica simulations, meaning
 that multiple instances (replicas) of your simulation are run
 simultaneously, with small amounts of data exchanged between replicas
 periodically.
 
 These are the relevant commands:
 
 "neb"_neb.html for nudged elastic band calculations
 "prd"_prd.html for parallel replica dynamics
 "temper"_temper.html for parallel tempering :ul
 
 NEB is a method for finding transition states and barrier energies.
 PRD is a method for performing accelerated dynamics to find and
 perform infrequent events.  Parallel tempering or replica exchange
 runs different replicas at a series of temperature to facilitate
 rare-event sampling.
 
 In all these cases, you must run with one or more processors per
 replica.  The processors assigned to each replica are determined at
 run-time by using the "-partition command-line
 switch"_Section_start.html#2_6 to launch LAMMPS on multiple
 partitions, which in this context are the same as replicas.  E.g.
 these commands:
 
 mpirun -np 16 lmp_linux -partition 8x2 -in in.temper
 mpirun -np 8 lmp_linux -partition 8x1 -in in.neb :pre
 
 would each run 8 replicas, on either 16 or 8 processors.  Note the use
 of the "-in command-line switch"_Section_start.html#2_6 to specify the
 input script which is required when running in multi-replica mode.
 
 Also note that with MPI installed on a machine (e.g. your desktop),
 you can run on more (virtual) processors than you have physical
 processors.  Thus the above commands could be run on a
 single-processor (or few-processor) desktop so that you can run
 a multi-replica simulation on more replicas than you have
 physical processors.
 
 :line
 
 4.6 Granular models :link(4_6),h4
 
 Granular system are composed of spherical particles with a diameter,
 as opposed to point particles.  This means they have an angular
 velocity and torque can be imparted to them to cause them to rotate.
 
 To run a simulation of a granular model, you will want to use
 the following commands:
 
 "atom_style"_atom_style.html granular
 "fix nve/sphere"_fix_nve_sphere.html
 "fix gravity"_fix_gravity.html :ul
 
 This compute
 
 "compute erotate/sphere"_compute_erotate_sphere.html :ul
 
 calculates rotational kinetic energy which can be "output with
 thermodynamic info"_Section_howto.html#4_15.
 
 Use one of these 3 pair potentials, which compute forces and torques
 between interacting pairs of particles:
 
 "pair_style"_pair_style.html gran/history
 "pair_style"_pair_style.html gran/no_history
 "pair_style"_pair_style.html gran/hertzian :ul
 
 These commands implement fix options specific to granular systems:
 
 "fix freeze"_fix_freeze.html
 "fix pour"_fix_pour.html
 "fix viscous"_fix_viscous.html
 "fix wall/gran"_fix_wall_gran.html :ul
 
 The fix style {freeze} zeroes both the force and torque of frozen
 atoms, and should be used for granular system instead of the fix style
 {setforce}.
 
 For computational efficiency, you can eliminate needless pairwise
 computations between frozen atoms by using this command:
 
 "neigh_modify"_neigh_modify.html exclude :ul
 
 :line
 
 4.7 TIP3P water model :link(4_7),h4
 
 The TIP3P water model as implemented in CHARMM
 "(MacKerell)"_#MacKerell specifies a 3-site rigid water molecule with
 charges and Lennard-Jones parameters assigned to each of the 3 atoms.
 In LAMMPS the "fix shake"_fix_shake.html command can be used to hold
 the two O-H bonds and the H-O-H angle rigid.  A bond style of
 {harmonic} and an angle style of {harmonic} or {charmm} should also be
 used.
 
 These are the additional parameters (in real units) to set for O and H
 atoms and the water molecule to run a rigid TIP3P-CHARMM model with a
 cutoff.  The K values can be used if a flexible TIP3P model (without
 fix shake) is desired.  If the LJ epsilon and sigma for HH and OH are
 set to 0.0, it corresponds to the original 1983 TIP3P model
 "(Jorgensen)"_#Jorgensen.
 
 O mass = 15.9994
 H mass = 1.008 :all(b),p
 
 O charge = -0.834
 H charge = 0.417 :all(b),p
 
 LJ epsilon of OO = 0.1521
 LJ sigma of OO = 3.1507
 LJ epsilon of HH = 0.0460
 LJ sigma of HH = 0.4000
 LJ epsilon of OH = 0.0836
 LJ sigma of OH = 1.7753 :all(b),p
 
 K of OH bond = 450
 r0 of OH bond = 0.9572 :all(b),p
 
 K of HOH angle = 55
 theta of HOH angle = 104.52 :all(b),p
 
 These are the parameters to use for TIP3P with a long-range Coulombic
 solver (Ewald or PPPM in LAMMPS), see "(Price)"_#Price for details:
 
 O mass = 15.9994
 H mass = 1.008 :all(b),p
 
 O charge = -0.830
 H charge = 0.415 :all(b),p
 
 LJ epsilon of OO = 0.102
 LJ sigma of OO = 3.188
 LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
 
 K of OH bond = 450
 r0 of OH bond = 0.9572 :all(b),p
 
 K of HOH angle = 55
 theta of HOH angle = 104.52 :all(b),p
 
 Wikipedia also has a nice article on "water
 models"_http://en.wikipedia.org/wiki/Water_model.
 
 :line
 
 4.8 TIP4P water model :link(4_8),h4
 
 The four-point TIP4P rigid water model extends the traditional
 three-point TIP3P model by adding an additional site, usually
 massless, where the charge associated with the oxygen atom is placed.
 This site M is located at a fixed distance away from the oxygen along
 the bisector of the HOH bond angle.  A bond style of {harmonic} and an
 angle style of {harmonic} or {charmm} should also be used.
 
 Currently, only a four-point model for long-range Coulombics is
 implemented via the LAMMPS "pair style
 lj/cut/coul/long/tip4p"_pair_lj.html.  A cutoff version may be added
 the future.  For both models, the bond lengths and bond angles should
 be held fixed using the "fix shake"_fix_shake.html command.
 
 These are the additional parameters (in real units) to set for O and H
 atoms and the water molecule to run a rigid TIP4P model with a cutoff
 "(Jorgensen)"_#Jorgensen.  Note that the OM distance is specified in
 the "pair_style"_pair_style.html command, not as part of the pair
 coefficients.
 
 O mass = 15.9994
 H mass = 1.008 :all(b),p
 
 O charge = -1.040
 H charge = 0.520 :all(b),p
 
 r0 of OH bond = 0.9572
 theta of HOH angle = 104.52 :all(b),p
 
 OM distance = 0.15 :all(b),p
 
 LJ epsilon of O-O = 0.1550
 LJ sigma of O-O = 3.1536
 LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
 
 These are the parameters to use for TIP4P with a long-range Coulombic
 solver (Ewald or PPPM in LAMMPS):
 
 O mass = 15.9994
 H mass = 1.008 :all(b),p
 
 O charge = -1.0484
 H charge = 0.5242 :all(b),p
 
 r0 of OH bond = 0.9572
 theta of HOH angle = 104.52 :all(b),p
 
 OM distance = 0.1250 :all(b),p
 
 LJ epsilon of O-O = 0.16275
 LJ sigma of O-O = 3.16435
 LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
 
 Wikipedia also has a nice article on "water
 models"_http://en.wikipedia.org/wiki/Water_model.
 
 :line
 
 4.9 SPC water model :link(4_9),h4
 
 The SPC water model specifies a 3-site rigid water molecule with
 charges and Lennard-Jones parameters assigned to each of the 3 atoms.
 In LAMMPS the "fix shake"_fix_shake.html command can be used to hold
 the two O-H bonds and the H-O-H angle rigid.  A bond style of
 {harmonic} and an angle style of {harmonic} or {charmm} should also be
 used.
 
 These are the additional parameters (in real units) to set for O and H
 atoms and the water molecule to run a rigid SPC model.
 
 O mass = 15.9994
 H mass = 1.008 :all(b),p
 
 O charge = -0.820
 H charge = 0.410 :all(b),p
 
 LJ epsilon of OO = 0.1553
 LJ sigma of OO = 3.166
 LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
 
 r0 of OH bond = 1.0
 theta of HOH angle = 109.47 :all(b),p
 
 Note that as originally proposed, the SPC model was run with a 9
 Angstrom cutoff for both LJ and Coulommbic terms.  It can also be used
 with long-range Coulombics (Ewald or PPPM in LAMMPS), without changing
 any of the parameters above, though it becomes a different model in
 that mode of usage.
 
 The SPC/E (extended) water model is the same, except
 the partial charge assignemnts change:
 
 O charge = -0.8476
 H charge = 0.4238 :all(b),p
 
 See the "(Berendsen)"_#Berendsen reference for more details on both
 the SPC and SPC/E models.
 
 Wikipedia also has a nice article on "water
 models"_http://en.wikipedia.org/wiki/Water_model.
 
 :line 
 
 4.10 Coupling LAMMPS to other codes :link(4_10),h4
 
 LAMMPS is designed to allow it to be coupled to other codes.  For
 example, a quantum mechanics code might compute forces on a subset of
 atoms and pass those forces to LAMMPS.  Or a continuum finite element
 (FE) simulation might use atom positions as boundary conditions on FE
 nodal points, compute a FE solution, and return interpolated forces on
 MD atoms.
 
 LAMMPS can be coupled to other codes in at least 3 ways.  Each has
 advantages and disadvantages, which you'll have to think about in the
 context of your application.
 
 (1) Define a new "fix"_fix.html command that calls the other code.  In
 this scenario, LAMMPS is the driver code.  During its timestepping,
 the fix is invoked, and can make library calls to the other code,
 which has been linked to LAMMPS as a library.  This is the way the
 "POEMS"_poems package that performs constrained rigid-body motion on
 groups of atoms is hooked to LAMMPS.  See the
 "fix_poems"_fix_poems.html command for more details.  See "this
 section"_Section_modify.html of the documentation for info on how to add
 a new fix to LAMMPS.
 
 :link(poems,http://www.rpi.edu/~anderk5/lab)
 
 (2) Define a new LAMMPS command that calls the other code.  This is
 conceptually similar to method (1), but in this case LAMMPS and the
 other code are on a more equal footing.  Note that now the other code
 is not called during the timestepping of a LAMMPS run, but between
 runs.  The LAMMPS input script can be used to alternate LAMMPS runs
 with calls to the other code, invoked via the new command.  The
 "run"_run.html command facilitates this with its {every} option, which
 makes it easy to run a few steps, invoke the command, run a few steps,
 invoke the command, etc.
 
 In this scenario, the other code can be called as a library, as in
 (1), or it could be a stand-alone code, invoked by a system() call
 made by the command (assuming your parallel machine allows one or more
 processors to start up another program).  In the latter case the
 stand-alone code could communicate with LAMMPS thru files that the
 command writes and reads.
 
 See "this section"_Section_modify.html of the documentation for how to
 add a new command to LAMMPS.
 
 (3) Use LAMMPS as a library called by another code.  In this case the
 other code is the driver and calls LAMMPS as needed.  Or a wrapper
 code could link and call both LAMMPS and another code as libraries.
 Again, the "run"_run.html command has options that allow it to be
 invoked with minimal overhead (no setup or clean-up) if you wish to do
 multiple short runs, driven by another program.
 
 Examples of driver codes that call LAMMPS as a library are included in
 the "couple" directory of the LAMMPS distribution; see couple/README
 for more details:
 
 simple: simple driver programs in C++ and C which invoke LAMMPS as a
 library :ulb,l
 
 lammps_quest: coupling of LAMMPS and "Quest"_quest, to run classical
 MD with quantum forces calculated by a density functional code :l
 
 lammps_spparks: coupling of LAMMPS and "SPPARKS"_spparks, to couple
 a kinetic Monte Carlo model for grain growth using MD to calculate
 strain induced across grain boundaries :l,ule
 
 :link(quest,http://dft.sandia.gov/Quest)
 :link(spparks,http://www.sandia.gov/~sjplimp/spparks.html)
 
 "This section"_Section_start.html#2_4 of the documentation describes
 how to build LAMMPS as a library.  Once this is done, you can
 interface with LAMMPS either via C++, C, or Fortran (or any other
 language that supports a vanilla C-like interface, e.g. a scripting
 language).  For example, from C++ you could create one (or more)
 "instances" of LAMMPS, pass it an input script to process, or execute
 individual commands, all by invoking the correct class methods in
 LAMMPS.  From C or Fortran you can make function calls to do the same
 things.  Library.cpp and library.h contain such a C interface with the
 functions:
 
 void lammps_open(int, char **, MPI_Comm, void **);
 void lammps_close(void *);
 void lammps_file(void *, char *);
 char *lammps_command(void *, char *); :pre
 
 The functions contain C++ code you could write in a C++ application
 that was invoking LAMMPS directly.  Note that LAMMPS classes are
 defined within a LAMMPS namespace (LAMMPS_NS) if you use them
 from another C++ application.
 
 Two of the routines in library.cpp are of particular note.  The
 lammps_open() function initiates LAMMPS and takes an MPI communicator
 as an argument.  It returns a pointer to a LAMMPS "object".  As with
 C++, the lammps_open() function can be called multiple times, to
 create multiple instances of LAMMPS.
 
 LAMMPS will run on the set of processors in the communicator.  This
 means the calling code can run LAMMPS on all or a subset of
 processors.  For example, a wrapper script might decide to alternate
 between LAMMPS and another code, allowing them both to run on all the
 processors.  Or it might allocate half the processors to LAMMPS and
 half to the other code and run both codes simultaneously before
 syncing them up periodically.
 
 Library.cpp contains a lammps_command() function to which the caller
 passes a single LAMMPS command (a string).  Thus the calling code can
 read or generate a series of LAMMPS commands (e.g. an input script)
 one line at a time and pass it thru the library interface to setup a
 problem and then run it.
 
 A few other sample functions are included in library.cpp, but the key
 idea is that you can write any functions you wish to define an
 interface for how your code talks to LAMMPS and add them to
 library.cpp and library.h.  The routines you add can access any LAMMPS
 data.  The examples/couple directory has example C++ and C codes which
 show how a stand-alone code can link LAMMPS as a library, run LAMMPS
 on a subset of processors, grab data from LAMMPS, change it, and put
 it back into LAMMPS.
 
 :line 
 
 4.11 Visualizing LAMMPS snapshots :link(4_11),h4
 
 LAMMPS itself does not do visualization, but snapshots from LAMMPS
 simulations can be visualized (and analyzed) in a variety of ways.
 
 LAMMPS snapshots are created by the "dump"_dump.html command which can
 create files in several formats.  The native LAMMPS dump format is a
 text file (see "dump atom" or "dump custom") which can be visualized
 by the "xmovie"_Section_tools.html#xmovie program, included with the
 LAMMPS package.  This produces simple, fast 2d projections of 3d
 systems, and can be useful for rapid debugging of simulation geometry
 and atom trajectories.
 
 Several programs included with LAMMPS as auxiliary tools can convert
 native LAMMPS dump files to other formats.  See the
 "Section_tools"_Section_tools.html doc page for details.  The first is
 the "ch2lmp tool"_Section_tools.html#charmm, which contains a
 lammps2pdb Perl script which converts LAMMPS dump files into PDB
 files.  The second is the "lmp2arc tool"_Section_tools.html#arc which
 converts LAMMPS dump files into Accelrys' Insight MD program files.
 The third is the "lmp2cfg tool"_Section_tools.html#cfg which converts
 LAMMPS dump files into CFG files which can be read into the
 "AtomEye"_atomeye visualizer.
 
 A Python-based toolkit distributed by our group can read native LAMMPS
 dump files, including custom dump files with additional columns of
 user-specified atom information, and convert them to various formats
 or pipe them into visualization software directly.  See the "Pizza.py
 WWW site"_pizza for details.  Specifically, Pizza.py can convert
 LAMMPS dump files into PDB, XYZ, "Ensight"_ensight, and VTK formats.
 Pizza.py can pipe LAMMPS dump files directly into the Raster3d and
 RasMol visualization programs.  Pizza.py has tools that do interactive
 3d OpenGL visualization and one that creates SVG images of dump file
 snapshots.
 
 LAMMPS can create XYZ files directly (via "dump xyz") which is a
 simple text-based file format used by many visualization programs
 including "VMD"_vmd.
 
 LAMMPS can create DCD files directly (via "dump dcd") which can be
 read by "VMD"_vmd in conjunction with a CHARMM PSF file.  Using this
 form of output avoids the need to convert LAMMPS snapshots to PDB
 files.  See the "dump"_dump.html command for more information on DCD
 files.
 
 LAMMPS can create XTC files directly (via "dump xtc") which is GROMACS
 file format which can also be read by "VMD"_vmd for visualization.
 See the "dump"_dump.html command for more information on XTC files.
 
 :link(pizza,http://www.sandia.gov/~sjplimp/pizza.html)
 :link(vmd,http://www.ks.uiuc.edu/Research/vmd)
 :link(ensight,http://www.ensight.com)
 :link(atomeye,http://mt.seas.upenn.edu/Archive/Graphics/A)
 
 :line
 
 4.12 Triclinic (non-orthogonal) simulation boxes :link(4_12),h4
 
 By default, LAMMPS uses an orthogonal simulation box to encompass the
 particles.  The "boundary"_boundary.html command sets the boundary
 conditions of the box (periodic, non-periodic, etc).  The orthogonal
 box has its "origin" at (xlo,ylo,zlo) and is defined by 3 edge vectors
 starting from the origin given by A = (xhi-xlo,0,0); B =
 (0,yhi-ylo,0); C = (0,0,zhi-zlo).  The 6 parameters
 (xlo,xhi,ylo,yhi,zlo,zhi) are defined at the time the simluation box
 is created, e.g. by the "create_box"_create_box.html or
 "read_data"_read_data.html or "read_restart"_read_restart.html
 commands.  Additionally, LAMMPS defines box size parameters lx,ly,lz
 where lx = xhi-xlo, and similarly in the y and z dimensions.  The 6
 parameters, as well as lx,ly,lz, can be output via the "thermo_style
 custom"_thermo_style.html command.
 
 LAMMPS also allows simulations to be perfored in non-orthogonal
 simulation boxes shaped as a parallelepiped with triclinic symmetry.
 The parallelepiped has its "origin" at (xlo,ylo,zlo) and is defined by
 3 edge vectors starting from the origin given by A = (xhi-xlo,0,0); B
 = (xy,yhi-ylo,0); C = (xz,yz,zhi-zlo).  {Xy,xz,yz} can be 0.0 or
 positive or negative values and are called "tilt factors" because they
 are the amount of displacement applied to faces of an originally
 orthogonal box to transform it into the parallelepiped.  Note that in
 LAMMPS the triclinic simulation box edge vectors A,B,C cannot be
 arbitrary vectors.  As indicated, A must be aligned with the x axis, B
 must be in the xy plane, and C is arbitrary.  However, this is not a
 restriction since it is possible to rotate any set of 3 crystal basis
 vectors so that they meet this restriction.
 
 The 9 parameters (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) are defined at the
 time the simluation box is created.  This happens in one of 3 ways.
 If the "create_box"_create_box.html command is used with a region of
 style {prism}, then a triclinic box is setup.  See the
 "region"_region.html command for details.  If the
 "read_data"_read_data.html command is used to define the simulation
 box, and the header of the data file contains a line with the "xy xz
 yz" keyword, then a triclinic box is setup.  See the
 "read_data"_read_data.html command for details.  Finally, if the
 "read_restart"_read_restart.html command reads a restart file which
 was written from a simulation using a triclinic box, then a triclinic
 box will be setup for the restarted simulation.
 
 Note that you can define a triclinic box with all 3 tilt factors =
 0.0, so that it is initially orthogonal.  This is necessary if the box
 will become non-orthogonal, e.g. due to the "fix npt"_fix_nh.html or
 "fix deform"_fix_deform.html commands.  Alternatively, you can use the
 "change_box"_change_box.html command to convert a simulation box from
 orthogonal to triclinic and vice versa.
 
 As with orthogonal boxes, LAMMPS defines triclinic box size parameters
 lx,ly,lz where lx = xhi-xlo, and similarly in the y and z dimensions.
 The 9 parameters, as well as lx,ly,lz, can be output via the
 "thermo_style custom"_thermo_style.html command.
 
 To avoid extremely tilted boxes (which would be computationally
 inefficient), no tilt factor can skew the box more than half the
 distance of the parallel box length, which is the 1st dimension in the
 tilt factor (x for xz).  For example, if xlo = 2 and xhi = 12, then
 the x box length is 10 and the xy tilt factor must be between -5 and
 5.  Similarly, both xz and yz must be between -(xhi-xlo)/2 and
 +(yhi-ylo)/2.  Note that this is not a limitation, since if the
 maximum tilt factor is 5 (as in this example), then configurations
 with tilt = ..., -15, -5, 5, 15, 25, ... are geometrically all
 equivalent.
 
 Triclinic crystal structures are often defined using three lattice
 constants {a}, {b}, and {c}, and three angles {alpha}, {beta} and
 {gamma}.  Note that in this nomenclature, the a,b,c lattice constants
 are the scalar lengths of the 3 A,B,C edge vectors defined above.  The
 relationship between these 6 quantities (a,b,c,alpha,beta,gamma) and
 the LAMMPS box sizes (lx,ly,lz) = (xhi-xlo,yhi-ylo,zhi-zlo) and tilt
 factors (xy,xz,yz) is as follows:
 
 :c,image(Eqs/box.jpg) 
 
 As discussed on the "dump"_dump.html command doc page, when the BOX
 BOUNDS for a snapshot is written to a dump file for a triclinic box,
 an orthogonal bounding box which encloses the triclinic simulation box
 is output, along with the 3 tilt factors (xy, xz, yz) of the triclinic
 box, formatted as follows:
 
 ITEM: BOX BOUNDS xy xz yz
 xlo_bound xhi_bound xy
 ylo_bound yhi_bound xz
 zlo_bound zhi_bound yz :pre
 
 This bounding box is convenient for many visualization programs and is
 calculated from the 9 triclinic box parameters
 (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) as follows:
 
 xlo_bound = xlo + MIN(0.0,xy,xz,xy+xz)
 xhi_bound = xhi + MAX(0.0,xy,xz,xy+xz)
 ylo_bound = ylo + MIN(0.0,yz)
 yhi_bound = yhi + MAX(0.0,yz)
 zlo_bound = zlo
 zhi_bound = zhi :pre
 
 These formulas can be inverted if you need to convert the bounding box
 back into the triclinic box parameters, e.g. xlo = xlo_bound -
 MIN(0.0,xy,xz,xy+xz).
 
 One use of triclinic simulation boxes is to model solid-state crystals
 with triclinic symmetry.  The "lattice"_lattice.html command can be
 used with non-orthogonal basis vectors to define a lattice that will
 tile a triclinic simulation box via the
 "create_atoms"_create_atoms.html command.
 
 A second use is to run Parinello-Rahman dyanamics via the "fix
 npt"_fix_nh.html command, which will adjust the xy, xz, yz tilt
 factors to compensate for off-diagonal components of the pressure
 tensor.  The analalog for an "energy minimization"_minimize.html is
 the "fix box/relax"_fix_box_relax.html command.
 
 A third use is to shear a bulk solid to study the response of the
 material.  The "fix deform"_fix_deform.html command can be used for
 this purpose.  It allows dynamic control of the xy, xz, yz tilt
 factors as a simulation runs.  This is discussed in the next section
 on non-equilibrium MD (NEMD) simulations.
 
 :line
 
 4.13 NEMD simulations :link(4_13),h4
 
 Non-equilibrium molecular dynamics or NEMD simulations are typically
 used to measure a fluid's rheological properties such as viscosity.
 In LAMMPS, such simulations can be performed by first setting up a
 non-orthogonal simulation box (see the preceding Howto section).
 
 A shear strain can be applied to the simulation box at a desired
 strain rate by using the "fix deform"_fix_deform.html command.  The
 "fix nvt/sllod"_fix_nvt_sllod.html command can be used to thermostat
 the sheared fluid and integrate the SLLOD equations of motion for the
 system.  Fix nvt/sllod uses "compute
 temp/deform"_compute_temp_deform.html to compute a thermal temperature
 by subtracting out the streaming velocity of the shearing atoms.  The
 velocity profile or other properties of the fluid can be monitored via
 the "fix ave/spatial"_fix_ave_spatial.html command.
 
 As discussed in the previous section on non-orthogonal simulation
 boxes, the amount of tilt or skew that can be applied is limited by
 LAMMPS for computational efficiency to be 1/2 of the parallel box
 length.  However, "fix deform"_fix_deform.html can continuously strain
 a box by an arbitrary amount.  As discussed in the "fix
 deform"_fix_deform.html command, when the tilt value reaches a limit,
 the box is re-shaped to the opposite limit which is an equivalent
 tiling of periodic space.  The strain rate can then continue to change
 as before.  In a long NEMD simulation these box re-shaping events may
 occur many times.
 
 In a NEMD simulation, the "remap" option of "fix
 deform"_fix_deform.html should be set to "remap v", since that is what
 "fix nvt/sllod"_fix_nvt_sllod.html assumes to generate a velocity
 profile consistent with the applied shear strain rate.
 
 An alternative method for calculating viscosities is provided via the
 "fix viscosity"_fix_viscosity.html command.
 
 :line
 
 4.14 Extended spherical and aspherical particles :link(4_14),h4
 
 Typical MD models treat atoms or particles as point masses.
 Sometimes, however, it is desirable to have a model with finite-size
 particles such as spherioids or aspherical ellipsoids.  The difference
 is that such particles have a moment of inertia, rotational energy,
 and angular momentum.  Rotation is induced by torque from interactions
 with other particles.
 
 LAMMPS has several options for running simulations with these kinds of
 particles.  The following aspects are discussed in turn:
 
 atom styles
 pair potentials
 time integration
 computes, thermodynamics, and dump output
 rigid bodies composed of extended particles :ul
 
 Atom styles :h5
 
 There are 3 "atom styles"_atom_style.html that allow for definition of
 finite-size particles: granular, dipole, ellipsoid.
 
 Granular particles are spheriods and each particle can have a unique
 diameter and mass (or density).  These particles store an angular
 velocity (omega) and can be acted upon by torque.
 
 Dipolar particles are typically spheriods with a point dipole and each
 particle type has a diamater and mass, set by the "shape"_shape.html
 and "mass"_mass.html commands.  These particles store an angular
 velocity (omega) and can be acted upon by torque.  They also store an
 orientation for the point dipole (mu) which has a length set by the
 "dipole"_dipole.html command.  The "set"_set.html command can be used
 to initialize the orientation of dipole moments.
 
 Ellipsoid particles are aspherical.  Each particle type has an
 ellipsoidal shape and mass, defined by the "shape"_shape.html and
 "mass"_mass.html commands.  These particles store an angular momentum
 and their orientation (quaternion), and can be acted upon by torque.
 They do not store an angular velocity (omega), which can be in a
 different direction than angular momentum, rather they compute it as
 needed.  Ellipsoidal particles can also store a dipole moment if an
 "atom_style hybrid ellipsoid dipole"_atom_style.html is used.  The
 "set"_set.html command can be used to initialize the orientation of
 ellipsoidal particles and has a brief explanation of quaternions.
 
 Note that if one of these atom styles is used (or multiple styles via
 the "atom_style hybrid"_atom_style.html command), not all particles in
 the system are required to be finite-size or aspherical.  For example,
 if the 3 shape parameters are set to the same value, the particle will
 be a spheroid rather than an ellipsoid.  If the 3 shape parameters are
 all set to 0.0 or if the diameter is set to 0.0, it will be a point
 particle.  If the dipole moment is set to zero, the particle will not
 have a point dipole associated with it.  The pair styles used to
 compute pairwise interactions will typically compute the correct
 interaction in these simplified (cheaper) cases.  "Pair_style
 hybrid"_pair_hybrid.html can be used to insure the correct
 interactions are computed for the appropriate style of interactions.
 Likewise, using groups to partition particles (ellipsoid versus
 spheroid versus point particles) will allow you to use the appropriate
 time integrators and temperature computations for each class of
 particles.  See the doc pages for various commands for details.
 
 Also note that for "2d simulations"_dimension.html, finite-size
 spheroids and ellipsoids are still treated as 3d particles, rather
 than as disks or ellipses.  This means they have the same moment of
 inertia for a 3d extended object.  When their temperature is
 coomputed, the correct degrees of freedom are used for rotation in a
 2d versus 3d system.
 
 Pair potentials :h5
 
 When a system with extended particles is defined, the particles will
 only rotate and experience torque if the force field computes such
 interactions.  These are the various "pair styles"_pair_style.html
 that generate torque:
 
 "pair_style gran/history"_pair_gran.html
 "pair_style gran/hertzian"_pair_gran.html
 "pair_style gran/no_history"_pair_gran.html
 "pair_style dipole/cut"_pair_dipole.html
 "pair_style gayberne"_pair_gayberne.html
 "pair_style resquared"_pair_resquared.html
 "pair_style lubricate"_pair_lubricate.html :ul
 
 The "granular pair styles"_pair_gran.html are used with "atom_style
 granular"_atom_style.html.  The "dipole pair style"_pair_dipole.html
 is used with "atom_style dipole"_atom_style.html.  The
 "GayBerne"_pair_gayberne.html and "REsquared"_pair_resquared.html
 potentials require particles have a "shape"_shape.html and are
 designed for "ellipsoidal particles"_atom_style.html.  The
 "lubrication potential"_pair_lubricate.html requires that particles
 have a "shape"_shape.html.  It can currently only be used with
 extended spherical particles.
 
 Time integration :h5
 
 There are 3 fixes that perform time integration on extended spherical
 particles, meaning the integrators update the rotational orientation
 and angular velocity or angular momentum of the particles:
 
 "fix nve/sphere"_fix_nve_sphere.html
 "fix nvt/sphere"_fix_nvt_sphere.html
 "fix npt/sphere"_fix_npt_sphere.html :ul
 
 Likewise, there are 3 fixes that perform time integration on extended
 aspherical particles:
 
 "fix nve/asphere"_fix_nve_asphere.html
 "fix nvt/asphere"_fix_nvt_asphere.html
 "fix npt/asphere"_fix_npt_asphere.html :ul
 
 The advantage of these fixes is that those which thermostat the
 particles include the rotational degrees of freedom in the temperature
 calculation and thermostatting.  Other thermostats can be used with
 fix nve/sphere or fix nve/asphere, such as fix langevin or fix
 temp/berendsen, but those thermostats only operate on the
 translational kinetic energy of the extended particles.
 
 Note that for mixtures of point and extended particles, you should
 only use these integration fixes on "groups"_group.html which contain
 extended particles.
 
 Computes, thermodynamics, and dump output :h5
 
 There are 4 computes that calculate the temperature or rotational energy
 of extended spherical or aspherical particles:
 
 "compute temp/sphere"_compute_temp_sphere.html
 "compute temp/asphere"_compute_temp_asphere.html
 "compute erotate/sphere"_compute_erotate_sphere.html
 "compute erotate/asphere"_compute_erotate_asphere.html :ul
 
 These include rotational degrees of freedom in their computation.  If
 you wish the thermodynamic output of temperature or pressure to use
 one of these computes (e.g. for a system entirely composed of extended
 particles), then the compute can be defined and the
 "thermo_modify"_thermo_modify.html command used.  Note that by
 default thermodynamic quantities will be calculated with a temperature
 that only includes translational degrees of freedom.  See the
 "thermo_style"_thermo_style.html command for details.
 
 The "dump custom"_dump.html command can output various attributes of
 extended particles, including the dipole moment (mu), the angular
 velocity (omega), the angular momentum (angmom), the quaternion
 (quat), and the torque (tq) on the particle.
 
 Rigid bodies composed of extended particles :h5
 
 The "fix rigid"_fix_rigid.html command treats a collection of
 particles as a rigid body, computes its inertia tensor, sums the total
 force and torque on the rigid body each timestep due to forces on its
 constituent particles, and integrates the motion of the rigid body.
 
 (NOTE: the feature described in the following paragraph has not yet
 been released.  It will be soon.)
 
 If any of the constituent particles of a rigid body are extended
 particles (spheroids or ellipsoids), then their contribution to the
 inertia tensor of the body is different than if they were point
 particles.  This means the rotational dynamics of the rigid body will
 be different.  Thus a model of a dimer is different if the dimer
 consists of two point masses versus two extended sphereoids, even if
 the two particles have the same mass.  Extended particles that
 experience torque due to their interaction with other particles will
 also impart that torque to a rigid body they are part of.
 
 See the "fix rigid" command for example of complex rigid-body models
 it is possible to define in LAMMPS.
 
 Note that the "fix shake"_fix_shake.html command can also be used to
 treat 2, 3, or 4 particles as a rigid body, but it always assumes the
 particles are point masses.
 
 :line
 
 4.15 Output from LAMMPS (thermo, dumps, computes, fixes, variables) :link(4_15),h4
 
 There are four basic kinds of LAMMPS output:
 
 "Thermodynamic output"_thermo_style.html, which is a list
 of quantities printed every few timesteps to the screen and logfile. :ulb,l
 
 "Dump files"_dump.html, which contain snapshots of atoms and various
 per-atom values and are written at a specified frequency. :l
 
 Certain fixes can output user-specified quantities to files: "fix
 ave/time"_fix_ave_time.html for time averaging, "fix
 ave/spatial"_fix_ave_spatial.html for spatial averaging, and "fix
 print"_fix_print.html for single-line output of
 "variables"_variable.html.  Fix print can also output to the
 screen. :l
 
 "Restart files"_restart.html. :l,ule
 
 A simulation prints one set of thermodynamic output and (optionally)
 restart files.  It can generate any number of dump files and fix
 output files, depending on what "dump"_dump.html and "fix"_fix.html
 commands you specify.
 
 As discussed below, LAMMPS gives you a variety of ways to determine
 what quantities are computed and printed when the thermodynamics,
 dump, or fix commands listed above perform output.  Throughout this
 discussion, note that users can also "add their own computes and fixes
 to LAMMPS"_Section_modify.html which can then generate values that can
 then be output with these commands.
 
 The following sub-sections discuss different LAMMPS command related
 to output and the kind of data they operate on and produce:
 
 "Global/per-atom/local data"_#global
 "Scalar/vector/array data"_#scalar
 "Thermodynamic output"_#thermo
 "Dump file output"_#dump
 "Fixes that write output files"_#fixoutput
 "Computes that process output quantities"_#computeoutput
 "Fixes that process output quantities"_#fixoutput
 "Computes that generate values to output"_#compute
 "Fixes that generate values to output"_#fix
 "Variables that generate values to output"_#variable
 "Summary table of output options and data flow between commands"_#table :ul
 
 Global/per-atom/local data :h5,link(global)
 
 Various output-related commands work with three different styles of
 data: global, per-atom, or local.  A global datum is one or more
 system-wide values, e.g. the temperature of the system.  A per-atom
 datum is one or more values per atom, e.g. the kinetic energy of each
 atom.  Local datums are calculated by each processor based on the
 atoms it owns, but there may be zero or more per atom, e.g. a list of
 bond distances.
 
 Scalar/vector/array data :h5,link(scalar)
 
 Global, per-atom, and local datums can each come in three kinds: a
 single scalar value, a vector of values, or a 2d array of values.  The
 doc page for a "compute" or "fix" or "variable" that generates data
 will specify both the style and kind of data it produces, e.g. a
 per-atom vector.
 
 When a quantity is accessed, as in many of the output commands
 discussed below, it can be referenced via the following bracket
 notation, where ID in this case is the ID of a compute.  The leading
 "c_" would be replaced by "f_" for a fix, or "v_" for a variable:
 
 c_ID | entire scalar, vector, or array
 c_ID\[I\] | one element of vector, one column of array
 c_ID\[I\]\[J\] | one element of array :tb(s=|)
 
 In other words, using one bracket reduces the dimension of the data
 once (vector -> scalar, array -> vector).  Using two brackets reduces
 the dimension twice (array -> scalar).  Thus a command that uses
 scalar values as input can typically also process elements of a vector
 or array.
 
 Thermodynamic output :h5,link(thermo)
 
 The frequency and format of thermodynamic output is set by the
 "thermo"_thermo.html, "thermo_style"_thermo_style.html, and
 "thermo_modify"_thermo_modify.html commands.  The
 "thermo_style"_thermo_style.html command also specifies what values
 are calculated and written out.  Pre-defined keywords can be specified
 (e.g. press, etotal, etc).  Three additional kinds of keywords can
 also be specified (c_ID, f_ID, v_name), where a "compute"_compute.html
 or "fix"_fix.html or "variable"_variable.html provides the value to be
 output.  In each case, the compute, fix, or variable must generate
 global values for input to the "thermo_style custom"_dump.html
 command.
 
 Dump file output :h5,link(dump)
 
 Dump file output is specified by the "dump"_dump.html and
 "dump_modify"_dump_modify.html commands.  There are several
 pre-defined formats (dump atom, dump xtc, etc).
 
 There is also a "dump custom"_dump.html format where the user
 specifies what values are output with each atom.  Pre-defined atom
 attributes can be specified (id, x, fx, etc).  Three additional kinds
 of keywords can also be specified (c_ID, f_ID, v_name), where a
 "compute"_compute.html or "fix"_fix.html or "variable"_variable.html
 provides the values to be output.  In each case, the compute, fix, or
 variable must generate per-atom values for input to the "dump
 custom"_dump.html command.
 
 There is also a "dump local"_dump.html format where the user specifies
 what local values to output.  A pre-defined index keyword can be
 specified to enumuerate the local values.  Two additional kinds of
 keywords can also be specified (c_ID, f_ID), where a
 "compute"_compute.html or "fix"_fix.html or "variable"_variable.html
 provides the values to be output.  In each case, the compute or fix
 must generate local values for input to the "dump local"_dump.html
 command.
 
 Fixes that write output files :h5,link(fixoutput)
 
 Sevarl fixes take various quantities as input and can write output
 files: "fix ave/time"_fix_ave_time.html, "fix
 ave/spatial"_fix_ave_spatial.html, "fix ave/histo"_fix_ave_histo.html,
 "fix ave/correlate"_fix_ave_correlate.html, and "fix
 print"_fix_print.html.
 
 The "fix ave/time"_fix_ave_time.html command enables direct output to
 a file and/or time-averaging of global scalars or vectors.  The user
 specifies one or more quantities as input.  These can be global
 "compute"_compute.html values, global "fix"_fix.html values, or
 "variables"_variable.html of any style except the atom style which
 produces per-atom values.  Since a variable can refer to keywords used
 by the "thermo_style custom"_thermo_style.html command (like temp or
 press) and individual per-atom values, a wide variety of quantities
 can be time averaged and/or output in this way.  If the inputs are one
 or more scalar values, then the fix generate a global scalar or vector
 of output.  If the inputs are one or more vector values, then the fix
 generates a global vector or array of output.  The time-averaged
 output of this fix can also be used as input to other output commands.
 
 The "fix ave/spatial"_fix_ave_spatial.html command enables direct
 output to a file of spatial-averaged per-atom quantities like those
 output in dump files, within 1d layers of the simulation box.  The
 per-atom quantities can be atom density (mass or number) or atom
 attributes such as position, velocity, force.  They can also be
 per-atom quantities calculated by a "compute"_compute.html, by a
 "fix"_fix.html, or by an atom-style "variable"_variable.html.  The
 spatial-averaged output of this fix can also be used as input to other
 output commands.
 
 The "fix ave/histo"_fix_ave_histo.html command enables direct output
 to a file of histogrammed quantities, which can be global or per-atom
 or local quantities.  The histogram output of this fix can also be
 used as input to other output commands.
 
 The "fix ave/correlate"_fix_ave_histo.html command enables direct
 output to a file of time-correlated quantities, which can be global
 scalars.  The correlation matrix output of this fix can also be used
 as input to other output commands.
 
 The "fix print"_fix_print.html command can generate a line of output
 written to the screen and log file or to a separate file, periodically
 during a running simulation.  The line can contain one or more
 "variable"_variable.html values for any style variable except the atom
 style).  As explained above, variables themselves can contain
 references to global values generated by "thermodynamic
 keywords"_thermo_style.html, "computes"_compute.html,
 "fixes"_fix.html, or other "variables"_variable.html, or to per-atom
 values for a specific atom.  Thus the "fix print"_fix_print.html
 command is a means to output a wide variety of quantities separate
 from normal thermodynamic or dump file output.
 
 Computes that process output quantities :h5,link(computeoutput)
 
 The "compute reduce"_compute_reduce.html and "compute
 reduce/region"_compute_reduce.html commands take one or more per-atom
 or local vector quantities as inputs and "reduce" them (sum, min, max,
 ave) to scalar quantities.  These are produced as output values which
 can be used as input to other output commands.
 
 The "compute property/atom"_compute_property_atom.html command takes a
 list of one or more pre-defined atom attributes (id, x, fx, etc) and
 stores the values in a per-atom vector or array.  These are produced
 as output values which can be used as input to other output commands.
 The list of atom attributes is the same as for the "dump
 custom"_dump.html command.
 
 The "compute property/local"_compute_property_local.html command takes
 a list of one or more pre-defined local attributes (bond info, angle
 info, etc) and stores the values in a local vector or array.  These
 are produced as output values which can be used as input to other
 output commands.
 
 The "compute atom/molecule"_compute_atom_molecule.html command takes a
 list of one or more per-atom quantities (from a compute, fix, per-atom
 variable) and sums the quantities on a per-molecule basis.  It
 produces a global vector or array as output values which can be used
 as input to other output commands.
 
 Fixes that process output quantities :h5,link(fixoutput)
 
 The "fix ave/atom"_fix_ave_atom.html command performs time-averaging
 of per-atom vectors.  The per-atom quantities can be atom attributes
 such as position, velocity, force.  They can also be per-atom
 quantities calculated by a "compute"_compute.html, by a
 "fix"_fix.html, or by an atom-style "variable"_variable.html.  The
 time-averaged per-atom output of this fix can be used as input to
 other output commands.
 
 The "fix store/state"_fix_store_state.html command can archive one or
 more per-atom attributes at a particular time, so that the old values
 can be used in a future calculation or output.  The list of atom
 attributes is the same as for the "dump custom"_dump.html command,
 including per-atom quantities calculated by a "compute"_compute.html,
 by a "fix"_fix.html, or by an atom-style "variable"_variable.html.
 The output of this fix can be used as input to other output commands.
 
 Computes that generate values to output :h5,link(compute)
 
 Every "compute"_compute.html in LAMMPS produces either global or
 per-atom or local values.  The values can be scalars or vectors or
 arrays of data.  These values can be output using the other commands
 described in this section.  The doc page for each compute command
 describes what it produces.  Computes that produce per-atom or local
 values have the word "atom" or "local" in their style name.  Computes
 without the word "atom" or "local" produce global values.
 
 Fixes that generate values to output :h5,link(fix)
 
 Some "fixes"_fix.html in LAMMPS produces either global or per-atom or
 local values which can be accessed by other commands.  The values can
 be scalars or vectors or arrays of data.  These values can be output
 using the other commands described in this section.  The doc page for
 each fix command tells whether it produces any output quantities and
 describes them.
 
 Variables that generate values to output :h5,link(variable)
 
 Every "variables"_variable.html defined in an input script generates
 either a global scalar value or a per-atom vector (only atom-style
 variables) when it is accessed.  The formulas used to define equal-
 and atom-style variables can contain references to the thermodynamic
 keywords and to global and per-atom data generated by computes, fixes,
 and other variables.  The values generated by variables can be output
 using the other commands described in this section.
 
 Summary table of output options and data flow between commands :h5,link(table)
 
 This table summarizes the various commands that can be used for
 generating output from LAMMPS.  Each command produces output data of
 some kind and/or writes data to a file.  Most of the commands can take
 data from other commands as input.  Thus you can link many of these
 commands together in pipeline form, where data produced by one command
 is used as input to another command and eventually written to the
 screen or to a file.  Note that to hook two commands together the
 output and input data types must match, e.g. global/per-atom/local
 data and scalar/vector/array data.
 
 Also note that, as described above, when a command takes a scalar as
 input, that could be an element of a vector or array.  Likewise a
 vector input could be a column of an array.
 
 Command: Input: Output:
 "thermo_style custom"_thermo_style.html: global scalars: screen, log file:
 "dump custom"_dump.html: per-atom vectors: dump file:
 "dump local"_dump.html: local vectors: dump file:
 "fix print"_fix_print.html: global scalar from variable: screen, file:
 "print"_print.html: global scalar from variable: screen:
 "computes"_compute.html: N/A: global/per-atom/local scalar/vector/array:
 "fixes"_fix.html: N/A: global/per-atom/local scalar/vector/array:
 "variables"_variable.html: global scalars, per-atom vectors: global scalar, per-atom vector:
 "compute reduce"_compute_reduce.html: per-atom/local vectors: global scalar/vector:
 "compute property/atom"_compute_property_atom.html: per-atom vectors: per-atom vector/array:
 "compute property/local"_compute_property_local.html: local vectors: local vector/array:
 "compute atom/molecule"_compute_atom_molecule.html: per-atom vectors: global vector/array:
 "fix ave/atom"_fix_ave_atom.html: per-atom vectors: per-atom vector/array:
 "fix ave/time"_fix_ave_time.html: global scalars/vectors: global scalar/vector/array, file:
 "fix ave/spatial"_fix_ave_spatial.html: per-atom vectors: global array, file:
 "fix ave/histo"_fix_ave_histo.html: global/per-atom/local scalars and vectors: global array, file:
 "fix ave/correlate"_fix_ave_correlate.html: global scalars: global array, file:
 "fix store/state"_fix_store_state.html: per-atom vectors: per-atom vector/array:
 :tb(s=:)
 
 :line
 
 4.16 Thermostatting, barostatting, and computing temperature :link(4_16),h4
 
 Thermostatting means controlling the temperature of particles in an MD
 simulation.  Barostatting means controlling the pressure.  Since the
 pressure includes a kinetic component due to particle velocities, both
 these operations require calculation of the temperature.  Typically a
 target temperature (T) and/or pressure (P) is specified by the user,
 and the thermostat or barostat attempts to equilibrate the system to
 the requested T and/or P.
 
 Temperature is computed as kinetic energy divided by some number of
 degrees of freedom (and the Boltzmann constant).  Since kinetic energy
 is a function of particle velocity, there is often a need to
 distinguish between a particle's advection velocity (due to some
 aggregate motiion of particles) and its thermal velocity.  The sum of
 the two is the particle's total velocity, but the latter is often what
 is wanted to compute a temperature.
 
 LAMMPS has several options for computing temperatures, any of which
 can be used in thermostatting and barostatting.  These "compute
 commands"_compute.html calculate temperature, and the "compute
 pressure"_compute_pressure.html command calculates pressure.
 
 "compute temp"_compute_temp.html
 "compute temp/sphere"_compute_temp_sphere.html
 "compute temp/asphere"_compute_temp_asphere.html
 "compute temp/com"_compute_temp_com.html
 "compute temp/deform"_compute_temp_deform.html
 "compute temp/partial"_compute_temp_partial.html
 "compute temp/profile"_compute_temp_profile.html
 "compute temp/ramp"_compute_temp_ramp.html
 "compute temp/region"_compute_temp_region.html :ul
 
 All but the first 3 calculate velocity biases (i.e. advection
 velocities) that are removed when computing the thermal temperature.
 "Compute temp/sphere"_compute_temp_sphere.html and "compute
 temp/asphere"_compute_temp_asphere.html compute kinetic energy for
 extended particles that includes rotational degrees of freedom.  They
 both allow, as an extra argument, which is another temperature compute
 that subtracts a velocity bias.  This allows the translational
 velocity of extended spherical or aspherical particles to be adjusted
 in prescribed ways.
 
 Thermostatting in LAMMPS is performed by "fixes"_fix.html, or in one
 case by a pair style.  Four thermostatting fixes are currently
 available: Nose-Hoover (nvt), Berendsen, Langevin, and direct
 rescaling (temp/rescale).  Dissipative particle dynamics (DPD)
 thermostatting can be invoked via the {dpd/tstat} pair style:
 
 "fix nvt"_fix_nh.html
 "fix nvt/sphere"_fix_nvt_sphere.html
 "fix nvt/asphere"_fix_nvt_asphere.html
 "fix nvt/sllod"_fix_nvt_sllod.html
 "fix temp/berendsen"_fix_temp_berendsen.html
 "fix langevin"_fix_langevin.html
 "fix temp/rescale"_fix_temp_rescale.html
 "pair_style dpd/tstat"_pair_dpd.html :ul
 
 "Fix nvt"_fix_nh.html only thermostats the translational velocity of
 particles.  "Fix nvt/sllod"_fix_nvt_sllod.html also does this, except
 that it subtracts out a velocity bias due to a deforming box and
 integrates the SLLOD equations of motion.  See the "NEMD
 simulations"_#4_13 section of this page for further details.  "Fix
 nvt/sphere"_fix_nvt_sphere.html and "fix
 nvt/asphere"_fix_nvt_asphere.html thermostat not only translation
 velocities but also rotational velocities for spherical and aspherical
 particles.
 
 DPD thermostatting alters pairwise interactions in a manner analagous
 to the per-particle thermostatting of "fix
 langevin"_fix_langevin.html.
 
 Any of the thermostatting fixes can use temperature computes that
 remove bias for two purposes: (a) computing the current temperature to
 compare to the requested target temperature, and (b) adjusting only
 the thermal temperature component of the particle's velocities.  See
 the doc pages for the individual fixes and for the
 "fix_modify"_fix_modify.html command for instructions on how to assign
 a temperature compute to a thermostatting fix.  For example, you can
 apply a thermostat to only the x and z components of velocity by using
 it in conjunction with "compute
 temp/partial"_compute_temp_partial.html.
 
 IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning
 they update the velocities and positions of particles due to forces
 and velocities respectively.  The other thermostat fixes only adjust
 velocities; they do NOT perform time integration updates.  Thus they
 should be used in conjunction with a constant NVE integration fix such
 as these:
 
 "fix nve"_fix_nve.html
 "fix nve/sphere"_fix_nve_sphere.html
 "fix nve/asphere"_fix_nve_asphere.html :ul
 
 Barostatting in LAMMPS is also performed by "fixes"_fix.html.  Two
 barosttating methods are currently available: Nose-Hoover (npt and
 nph) and Berendsen:
 
 "fix npt"_fix_nh.html
 "fix npt/sphere"_fix_npt_sphere.html
 "fix npt/asphere"_fix_npt_asphere.html
 "fix nph"_fix_nh.html
 "fix press/berendsen"_fix_press_berendsen.html :ul
 
 The "fix npt"_fix_nh.html commands include a Nose-Hoover thermostat
 and barostat.  "Fix nph"_fix_nh.html is just a Nose/Hoover barostat;
 it does no thermostatting.  Both "fix nph"_fix_nh.html and "fix
 press/bernendsen"_fix_press_berendsen.html can be used in conjunction
 with any of the thermostatting fixes.
 
 As with the thermostats, "fix npt"_fix_nh.html and "fix
 nph"_fix_nh.html only use translational motion of the particles in
 computing T and P and performing thermo/barostatting.  "Fix
 npt/sphere"_fix_npt_sphere.html and "fix
 npt/asphere"_fix_npt_asphere.html thermo/barostat using not only
 translation velocities but also rotational velocities for spherical
 and aspherical particles.
 
 All of the barostatting fixes use the "compute
 pressure"_compute_pressure.html compute to calculate a current
 pressure.  By default, this compute is created with a simple "compute
 temp"_compute_temp.html (see the last argument of the "compute
 pressure"_compute_pressure.html command), which is used to calculated
 the kinetic componenet of the pressure.  The barostatting fixes can
 also use temperature computes that remove bias for the purpose of
 computing the kinetic componenet which contributes to the current
 pressure.  See the doc pages for the individual fixes and for the
 "fix_modify"_fix_modify.html command for instructions on how to assign
 a temperature or pressure compute to a barostatting fix.
 
 IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods ("fix
 npt"_fix_nh.html and "fix nph"_fix_nh.html) perform time
 integration.  "Fix press/berendsen"_fix_press_berendsen.html does NOT,
 so it should be used with one of the constant NVE fixes or with one of
 the NVT fixes.
 
 Finally, thermodynamic output, which can be setup via the
 "thermo_style"_thermo_style.html command, often includes temperature
 and pressure values.  As explained on the doc page for the
 "thermo_style"_thermo_style.html command, the default T and P are
 setup by the thermo command itself.  They are NOT the ones associated
 with any thermostatting or barostatting fix you have defined or with
 any compute that calculates a temperature or pressure.  Thus if you
 want to view these values of T and P, you need to specify them
 explicitly via a "thermo_style custom"_thermo_style.html command.  Or
 you can use the "thermo_modify"_thermo_modify.html command to
 re-define what temperature or pressure compute is used for default
 thermodynamic output.
 
 :line
 
 4.17 Walls :link(4_17),h4
 
 Walls in an MD simulation are typically used to bound particle motion,
 i.e. to serve as a boundary condition.
 
 Walls in LAMMPS can be of rough (made of particles) or idealized
 surfaces.  Ideal walls can be smooth, generating forces only in the
 normal direction, or frictional, generating forces also in the
 tangential direction.
 
 Rough walls, built of particles, can be created in various ways.  The
 particles themselves can be generated like any other particle, via the
 "lattice"_lattice.html and "create_atoms"_create_atoms.html commands,
 or read in via the "read_data"_read_data.html command.
 
 Their motion can be constrained by many different commands, so that
 they do not move at all, move together as a group at constant velocity
 or in response to a net force acting on them, move in a prescribed
 fashion (e.g. rotate around a point), etc.  Note that if a time
 integration fix like "fix nve"_fix_nve.html or "fix nvt"_fix_nh.html
 is not used with the group that contains wall particles, their
 positions and velocities will not be updated.
 
 "fix aveforce"_fix_aveforce.html - set force on particles to average value, so they move together
 "fix setforce"_fix_setforce.html - set force on particles to a value, e.g. 0.0
 "fix freeze"_fix_freeze.html - freeze particles for use as granular walls
 "fix nve/noforce"_fix_nve_noforce.html - advect particles by their velocity, but without force
 "fix move"_fix_move.html - prescribe motion of particles by a linear velocity, oscillation, rotation, variable :ul
 
 The "fix move"_fix_move.html command offers the most generality, since
 the motion of individual particles can be specified with
 "variable"_variable.html formula which depends on time and/or the
 particle position.
 
 For rough walls, it may be useful to turn off pairwise interactions
 between wall particles via the "neigh_modify
 exclude"_neigh_modify.html command.
 
 Rough walls can also be created by specifying frozen particles that do
 not move and do not interact with mobile particles, and then tethering
 other particles to the fixed particles, via a "bond"_bond_style.html.
 The bonded particles do interact with other mobile particles.
 
 Idealized walls can be specified via several fix commands.  "Fix
 wall/gran"_fix_wall_gran.html creates frictional walls for use with
 granular particles; all the other commands create smooth walls.
 
 "fix wall/reflect"_fix_wall_reflect.html - reflective flat walls
 "fix wall/lj93"_fix_wall.html - flat walls, with Lennard-Jones 9/3 potential
 "fix wall/lj126"_fix_wall.html - flat walls, with Lennard-Jones 12/6 potential
 "fix wall/colloid"_fix_wall.html - flat walls, with "pair_style colloid"_pair_colloid.html potential
 "fix wall/harmonic"_fix_wall.html - flat walls, with repulsive harmonic spring potential
 "fix wall/region"_fix_wall_region.html - use region surface as wall
 "fix wall/gran"_fix_wall_gran.html - flat or curved walls with "pair_style granular"_pair_gran.html potential :ul
 
 The {lj93}, {lj126}, {colloid}, and {harmonic} styles all allow the
 flat walls to move with a constant velocity, or oscillate in time.
 The "fix wall/region"_fix_wall_region.html command offers the most
 generality, since the region surface is treated as a wall, and the
 geometry of the region can be a simple primitive volume (e.g. a
 sphere, or cube, or plane), or a complex volume made from the union
 and intersection of primitive volumes.  "Regions"_region.html can also
 specify a volume "interior" or "exterior" to the specified primitive
 shape or {union} or {intersection}.  "Regions"_region.html can also be
 "dynamic" meaning they move with constant velocity, oscillate, or
 rotate.
 
 The only frictional idealized walls currently in LAMMPS are flat or
 curved surfaces specified by the "fix wall/gran"_fix_wall_gran.html
 command.  At some point we plan to allow regoin surfaces to be used as
 frictional walls, as well as triangulated surfaces.
 
 :line
 
 4.18 Elastic constants :link(4_18),h4
 
 Elastic constants characterize the stiffness of a material. The formal
 definition is provided by the linear relation that holds between the
 stress and strain tensors in the limit of infinitesimal deformation.
 In tensor notation, this is expressed as s_ij = C_ijkl * e_kl, where
 the repeated indices imply summation. s_ij are the elements of the
 symmetric stress tensor. e_kl are the elements of the symmetric strain
 tensor. C_ijkl are the elements of the fourth rank tensor of elastic
 constants. In three dimensions, this tensor has 3^4=81 elements. Using
 Voigt notation, the tensor can be written as a 6x6 matrix, where C_ij
 is now the derivative of s_i w.r.t. e_j. Because s_i is itself a
 derivative w.r.t. e_i, it follows that C_ij is also symmetric, with at
 most 7*6/2 = 21 distinct elements.
 
 At zero temperature, it is easy to estimate these derivatives by
 deforming the cell in one of the six directions using the command
 "displace_box"_displace_box.html and measuring the change in the
 stress tensor. A general-purpose script that does this is given in the
 examples/elastic directory described in "this
 section"_Section_example.html.
 
 Calculating elastic constants at finite temperature is more
 challenging, because it is necessary to run a simulation that perfoms
 time averages of differential properties. One way to do this is to
 measure the change in average stress tensor in an NVT simulations when
 the cell volume undergoes a finite deformation. In order to balance
 the systematic and statistical errors in this method, the magnitude of
 the deformation must be chosen judiciously, and care must be taken to
 fully equilibrate the deformed cell before sampling the stress
 tensor. Another approach is to sample the triclinic cell fluctuations
 that occur in an NPT simulation. This method can also be slow to
 converge and requires careful post-processing "(Shinoda)"_#Shinoda
 
-:line
-
-4.19 Computing free energies from thermodynamic integration :link(4_19),h4
-
-Thermodynamic integration is a widely used method to compute free
-energies from atomistic simulations.  LAMMPS can be used to run
-thermodynamic integration calculations using the methods discussed in
-this section and the "fix adapt"_fix_adapt.html command.  Currently,
-it is capable of the transformations essential for computing melting
-points using the pseudo-supercritical path method developed by "Eike
-and Maginn"_#Eike_Maginn.
-
-See the examples/TI directory for more information and sample files
-that compute a melting point using the techniques described in this
-section.  That directory has its own README file.  See also the paper
-by "Jayaraman"_#Jayaraman for an example of using this implementation
-of thermodynamic integration in LAMMPS to compute melting points of
-alkali nitrate salts, using the steps outlined here.
-
-In this method, three intermediate "pseudo-supercritical" states are
-accessed in the transformation between the liquid and solid
-phases. These pseudo-states are a weak liquid, a dense weak liquid,
-and an ordered weak phase. The transformation between the liquid and
-solid states can also be driven uisng the "fix adapt"_fix_adapt.html
-command.
-
-For the transformation from the liquid to the weak phase, the
-intermolecular interactions need to be weakened.  Appropriate scale
-factors, computed by variables you define, and applied to pair styles
-by "fix adapt"_fix_adapt.html, can be used to do this, as in the
-example scripts.  The "compute ti"_compute_ti.html command can
-accumulate the value of dU/d{lambda}.  See "Jayaraman and
-Maginn"_#Jayaraman_Maginn for more information about calculating a
-free energy from dU/d{lambda}.
- 
-IMPORTANT NOTE: The pair styles that fix adapt can scale on-the-fly
-are listed on the "fix adapt"_fix_adapt doc page.  interaction scaling
-is desired.  If a pair style is not on that list, it is generally
-quite easy to add an extract() method to the pair style, to enable fix
-adapt to rescale it.
-
-Step 2 is the transformation of the simulation box density from the
-liquid phase to that of the equilibrated crystal.  The parameters for
-box1 and box2 should be obtained from equilibrated NPT simulations of
-the liquid and crystal phases and used in a "fix
-deform"_fix_deform.html command to change the box size and/or shape.
-It also advisable to use "fix adapt"_fix_adapt.html on the pair styles
-to prevent overlaps which may occur during the box transformation.
-
-In step 3, the dense, weak system is transformed to an ordered state,
-which has the same ordering as in the equilibrated crystal.  Ordering
-is achieved by introducing an attractive potential between atoms and
-lattice sites.  These lattice sites can be calculated as the mean
-positions of the atoms in an equilibrium simulation of the
-crystal. The "pair/gauss"_pair_gauss.html command can be used to
-introduce an attractive Gaussian potential between the atoms and their
-corresponding lattice sites. The prefactor of the Gaussian pair
-potential can be scaled by "fix adapt"_fix_adapt.html to turn on the
-attractions.  Again, the quantity dU/d{lambda} can be tracked via the
-"compute ti"_compute_ti.html command.
-
-Step 4 is the transformation of the ordered state to the final
-crystal.  In this step, the intermolecular interactions are scaled
-back to full strength, while the Gaussian tethers are removed, all via
-"fix adapt"_fix_adapt.html.
-
 :line
 :line
 
 :link(Berendsen)
 [(Berendsen)] Berendsen, Grigera, Straatsma, J Phys Chem, 91,
 6269-6271 (1987).
 
 :link(Cornell)
 [(Cornell)] Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
 Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
 
 :link(Horn)
 [(Horn)] Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
 J Chem Phys, 120, 9665 (2004).
 
 :link(MacKerell)
 [(MacKerell)] MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
 Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
 
+:link(Mayo)
+[(Mayo)] Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
+(1990).
+
 :link(Jorgensen)
 [(Jorgensen)] Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
 Phys, 79, 926 (1983).
 
 :link(Price)
 [(Price)] Price and Brooks, J Chem Phys, 121, 10096 (2004).
 
 :link(Shinoda)
 [(Shinoda)] Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
-
-:link(Eike_Maginn) 
-[(Eike and Maginn)] Eike and Maginn, J Chem Phys, 124,
-164503 (2006). 
-
-:link(Jayaraman_Maginn)
-[(Jayaraman and Maginn)] Jayaraman and Maginn, Journal of Chemical Physics, 
-127, 214504 (2007).
-
-:link(Jayaraman)
-[(Jayaraman)] Jayaraman, Thompson, von Lilienfeld and Maginn, Industrial 
-and Engineering Chemistry Research, 49, 559-571 (2010).   
diff --git a/doc/compute_pair.html b/doc/compute_pair.html
index 0a00afd19..fbb946c8d 100644
--- a/doc/compute_pair.html
+++ b/doc/compute_pair.html
@@ -1,69 +1,85 @@
 <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>compute pair command 
 </H3>
 <P><B>Syntax:</B>
 </P>
-<PRE>compute ID group-ID pair pstyle 
+<PRE>compute ID group-ID pair pstyle evalue 
 </PRE>
 <UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
 <LI>pair = style name of this compute command
-<LI>pstyle = style name of a pair style that calculates additional values 
+<LI>pstyle = style name of a pair style that calculates additional values
+<LI>evalue = <I>epair</I> or <I>evdwl</I> or <I>evoul</I> or blank (optional setting) 
 </UL>
 <P><B>Examples:</B>
 </P>
 <PRE>compute 1 all pair gauss
+compute 1 all pair lj/cut/coul/cut ecoul
 compute 1 all pair reax 
 </PRE>
 <P><B>Description:</B>
 </P>
 <P>Define a computation that extracts additional values calculated by a
 pair style, sums them across processors, and makes them accessible for
 output or further processing by other commands.  The group specified
 for this command is ignored.
 </P>
-<P>The specified <I>pstyle</I> must be a pair style that produces additional
-values.  If a <A HREF = "pair_hybrid.html">hybrid pair style</A> is used, then
-<I>pstyle</I> should be the name of a sub-style.
+<P>The specified <I>pstyle</I> must be a pair style used in your simulation
+either by itself or as a sub-style in a <A HREF = "pair_hybrid.html">pair_style hybrid or
+hybrid/overlay</A> command.
 </P>
-<P>All pair styles tally a potential energy, which is accessed by the
-<A HREF = "compute_pe.html">compute pe</A> and <A HREF = "compute_pe_atom.html">compute
-pe/atom</A> commands.  Some pair styles tally one or
-more additional values, such as a breakdown of the total pair
-potential energy into sub-categories.  See the doc page for
-<A HREF = "pair_style.html">individual pair styles</A> for info on these values.
+<P>The <I>evalue</I> setting is optional; it may be left off the command.  All
+pair styles tally a potential energy <I>epair</I> which may be broken into
+two parts: <I>evdwl</I> and <I>ecoul</I> such that <I>epair</I> = <I>evdwl</I> + <I>evoul</I>.
+If the pair style calculates Coulombic interactions, their energy will
+be tallied in <I>ecoul</I>.  Everything else (whether it is a Lennard-Jones
+style van der Waals interaction or not) is tallied in <I>evdwl</I>.  If
+<I>evalue</I> is specified as <I>epair</I> or left out, then <I>epair</I> is stored
+as a global scalar by this compute.  This is useful when using
+<A HREF = "pair_hybrid.html">pair_style hybrid</A> if you want to know the portion
+of the total energy contributed by one sub-style.  If <I>evalue</I> is
+specfied as <I>evdwl</I> or <I>ecoul</I>, then just that portion of the energy
+is stored as a global scalar.
 </P>
-<P>The compute pair command lets you access this data as a global vector
-of values and then use other <A HREF = "Section_howto.html#4_15">output options</A>
-that work with <A HREF = "compute.html">compute commands</A> to see or use the
-values.
+<P>Some pair styles tally additional quantities, e.g. a breakdown of
+potential energy into a dozen or so components is tallied by the
+<A HREF = "pair_reax.html">pair_style reax</A> commmand.  These values (1 or more)
+are stored as a global vector by this compute.  See the doc page for
+<A HREF = "pair_style.html">individual pair styles</A> for info on these values.
 </P>
 <P><B>Output info:</B>
 </P>
-<P>This compute calculates a global vector of length >= 1, as determined
-by the pair style.  These values can be used by any command that uses
-global vector values from a compute as input.  See <A HREF = "Section_howto.html#4_15">this
-section</A> for an overview of LAMMPS output
-options.
+<P>This compute calculates a global scalar which is <I>epair</I> or <I>evdwl</I> or
+<I>evoul</I>.  If the pair style supports it, it also calculates a global
+vector of length >= 1, as determined by the pair style.  These values
+can be used by any command that uses global scalar or vector values
+from a compute as input.  See <A HREF = "Section_howto.html#4_15">this section</A>
+for an overview of LAMMPS output options.
 </P>
-<P>The vector values calculated by this compute are "extensive".  They
-are in whatever units the pair style produces.
+<P>The scalar and vector values calculated by this compute are
+"extensive".
+</P>
+<P>The scalar value will be in energy <A HREF = "units.html">units</A>.  The vector
+values will typically also be in energy <A HREF = "units.html">units</A>, but
+see the doc page for the pair style for details.
 </P>
 <P><B>Restrictions:</B> none
 </P>
 <P><B>Related commands:</B>
 </P>
 <P><A HREF = "compute_pe.html">compute pe</A>
 </P>
-<P><B>Default:</B> none
+<P><B>Default:</B>
+</P>
+<P>The default for <I>evalue</I> is <I>epair</I>.
 </P>
 </HTML>
diff --git a/doc/compute_pair.txt b/doc/compute_pair.txt
index 26b2e4dd8..3cca4656c 100644
--- a/doc/compute_pair.txt
+++ b/doc/compute_pair.txt
@@ -1,64 +1,80 @@
 "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
 
 compute pair command :h3
 
 [Syntax:]
 
-compute ID group-ID pair pstyle :pre
+compute ID group-ID pair pstyle evalue :pre
 
 ID, group-ID are documented in "compute"_compute.html command
 pair = style name of this compute command
-pstyle = style name of a pair style that calculates additional values :ul
+pstyle = style name of a pair style that calculates additional values
+evalue = {epair} or {evdwl} or {evoul} or blank (optional setting) :ul
 
 [Examples:]
 
 compute 1 all pair gauss
+compute 1 all pair lj/cut/coul/cut ecoul
 compute 1 all pair reax :pre
 
 [Description:]
 
 Define a computation that extracts additional values calculated by a
 pair style, sums them across processors, and makes them accessible for
 output or further processing by other commands.  The group specified
 for this command is ignored.
 
-The specified {pstyle} must be a pair style that produces additional
-values.  If a "hybrid pair style"_pair_hybrid.html is used, then
-{pstyle} should be the name of a sub-style.
-
-All pair styles tally a potential energy, which is accessed by the
-"compute pe"_compute_pe.html and "compute
-pe/atom"_compute_pe_atom.html commands.  Some pair styles tally one or
-more additional values, such as a breakdown of the total pair
-potential energy into sub-categories.  See the doc page for
+The specified {pstyle} must be a pair style used in your simulation
+either by itself or as a sub-style in a "pair_style hybrid or
+hybrid/overlay"_pair_hybrid.html command.
+
+The {evalue} setting is optional; it may be left off the command.  All
+pair styles tally a potential energy {epair} which may be broken into
+two parts: {evdwl} and {ecoul} such that {epair} = {evdwl} + {evoul}.
+If the pair style calculates Coulombic interactions, their energy will
+be tallied in {ecoul}.  Everything else (whether it is a Lennard-Jones
+style van der Waals interaction or not) is tallied in {evdwl}.  If
+{evalue} is specified as {epair} or left out, then {epair} is stored
+as a global scalar by this compute.  This is useful when using
+"pair_style hybrid"_pair_hybrid.html if you want to know the portion
+of the total energy contributed by one sub-style.  If {evalue} is
+specfied as {evdwl} or {ecoul}, then just that portion of the energy
+is stored as a global scalar.
+
+Some pair styles tally additional quantities, e.g. a breakdown of
+potential energy into a dozen or so components is tallied by the
+"pair_style reax"_pair_reax.html commmand.  These values (1 or more)
+are stored as a global vector by this compute.  See the doc page for
 "individual pair styles"_pair_style.html for info on these values.
 
-The compute pair command lets you access this data as a global vector
-of values and then use other "output options"_Section_howto.html#4_15
-that work with "compute commands"_compute.html to see or use the
-values.
-
 [Output info:]
 
-This compute calculates a global vector of length >= 1, as determined
-by the pair style.  These values can be used by any command that uses
-global vector values from a compute as input.  See "this
-section"_Section_howto.html#4_15 for an overview of LAMMPS output
-options.
+This compute calculates a global scalar which is {epair} or {evdwl} or
+{evoul}.  If the pair style supports it, it also calculates a global
+vector of length >= 1, as determined by the pair style.  These values
+can be used by any command that uses global scalar or vector values
+from a compute as input.  See "this section"_Section_howto.html#4_15
+for an overview of LAMMPS output options.
 
-The vector values calculated by this compute are "extensive".  They
-are in whatever units the pair style produces.
+The scalar and vector values calculated by this compute are
+"extensive".
+
+The scalar value will be in energy "units"_units.html.  The vector
+values will typically also be in energy "units"_units.html, but
+see the doc page for the pair style for details.
 
 [Restrictions:] none
 
 [Related commands:]
 
 "compute pe"_compute_pe.html
 
-[Default:] none
+[Default:]
+
+The default for {evalue} is {epair}.
diff --git a/doc/special_bonds.html b/doc/special_bonds.html
index 3a6326b1e..805f409c9 100644
--- a/doc/special_bonds.html
+++ b/doc/special_bonds.html
@@ -1,153 +1,182 @@
 <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>special_bonds command 
 </H3>
 <P><B>Syntax:</B>
 </P>
 <PRE>special_bonds keyword values ... 
 </PRE>
 <UL><LI>one or more keyword/value pairs may be appended 
 
-<LI>keyword = <I>amber</I> or <I>charmm</I> or <I>fene</I> or <I>lj/coul</I> or <I>lj</I> or <I>coul</I> or <I>angle</I> or <I>dihedral</I> or <I>extra</I> 
+<LI>keyword = <I>amber</I> or <I>charmm</I> or <I>dreiding</I> or <I>fene</I> or <I>lj/coul</I> or <I>lj</I> or <I>coul</I> or <I>angle</I> or <I>dihedral</I> or <I>extra</I> 
 
 <PRE>  <I>amber</I> values = none
   <I>charmm</I> values = none
+  <I>dreiding</I> values = none
   <I>fene</I> values = none
   <I>lj/coul</I> values = w1,w2,w3
     w1,w2,w3 = weights (0.0 to 1.0) on pairwise Lennard-Jones and Coulombic interactions
   <I>lj</I> values = w1,w2,w3
     w1,w2,w3 = weights (0.0 to 1.0) on pairwise Lennard-Jones interactions
   <I>coul</I> values = w1,w2,w3
     w1,w2,w3 = weights (0.0 to 1.0) on pairwise Coulombic interactions
   <I>angle</I> value = <I>yes</I> or <I>no</I>
   <I>dihedral</I> value = <I>yes</I> or <I>no</I>
   <I>extra</I> value = N
     N = number of extra 1-2,1-3,1-4 interactions to save space for 
 </PRE>
 
 </UL>
 <P>Examples:
 </P>
 <PRE>special_bonds amber
 special_bonds charmm
 special_bonds fene dihedral no
 special_bonds lj/coul 0.0 0.0 0.5 angle yes dihedral yes
 special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.0 dihedral yes
 special_bonds lj/coul 0 1 1 extra 2 
 </PRE>
 <P><B>Description:</B>
 </P>
 <P>Set weighting coefficients for pairwise energy and force contributions
 from atom pairs that are also bonded to each other directly or
 indirectly.  For Lennard-Jones (LJ) and Coulombic pairwise
 interactions, these coefficients come in sets of three.  The 1st
 coefficient is the weighting factor on 1-2 atom pairs, which are those
 directly bonded to each other.  The 2nd coefficient is the weighting
 factor on 1-3 atom pairs which are those separated by 2 bonds
 (e.g. the two H atoms in a water molecule).  The 3rd coefficient is
 the weighting factor on 1-4 atom pairs which are those separated by 3
 bonds (e.g. the 1st and 4th atoms in a dihedral interaction).  Thus if
 the 1-2 coefficient is set to 0.0, then the pairwise interaction is
 effectively turned off for all pairs of atoms bonded to each other.
 </P>
 <P>IMPORTANT NOTE: For purposes of computing weighted pairwise
 interactions, 1-3 and 1-4 interactions are not defined from the list
 of angles or dihedrals used by the simulation.  Rather, they are
 inferred topologically from the set of bonds defined when the
 simulation is defined from a data or restart file (see
 <A HREF = "read_data.html">read_data</A> or <A HREF = "read_restart.html">read_restart</A>
 commands).  Thus the set of 1-2,1-3,1-4 interactions that the weights
 apply to is the same whether angle and dihedral potentials are
 computed or not, and remains the same even if bonds are constrained,
 or turned off, or removed during a simulation.
 </P>
 <P>The two exceptions to this rule are (a) if the <I>angle</I> or <I>dihedral</I>
 keywords are set to <I>yes</I> (see below), or (b) if the
 <A HREF = "delete_bonds.html">delete_bonds</A> command is used with the <I>special</I>
 option that recomputes the 1-2,1-3,1-4 topologies after bonds are
 deleted; see the <A HREF = "delete_bonds.html">delete_bonds</A> command for more
 details.
 </P>
 <P>The <I>amber</I> keyword sets the 3 coefficients to 0.0, 0.0, 0.5 for LJ
 interactions and to 0.0, 0.0, 0.8333 for Coulombic interactions, which
 is the default for a commonly used version of the AMBER force field,
-where the last value is really 5/6.
+where the last value is really 5/6.  See <A HREF = "#Cornell">(Cornell)</A> for a
+description of the AMBER force field.
 </P>
 <P>The <I>charmm</I> keyword sets the 3 coefficients to 0.0, 0.0, 0.0 for both
 LJ and Coulombic interactions, which is the default for a commonly
 used version of the CHARMM force field.  Note that in pair styles
 <I>lj/charmm/coul/charmm</I> and <I>lj/charmm/coul/long</I> the 1-4 coefficients
 are defined explicitly, and these pairwise contributions are computed
 as part of the charmm dihedral style - see the
 <A HREF = "pair_coeff.html">pair_coeff</A> and <A HREF = "dihedral_style.html">dihedral_style</A>
-commands for more information.
+commands for more information.  See <A HREF = "#MacKerell">(MacKerell)</A> for a
+description of the CHARMM force field.
+</P>
+<P>The <I>dreiding</I> keyword sets the 3 coefficients to 0.0, 0.0, 1.0 for both
+LJ and Coulombic interactions, which is the default for the Dreiding
+force field, as discussed in <A HREF = "#Mayo">(Mayo)</A>.
 </P>
 <P>The <I>fene</I> keyword sets the 3 coefficients to 0.0, 1.0, 1.0 for both
 LJ and Coulombic interactions, which is consistent with a
-coarse-grained polymer model with <A HREF = "bond_fene.html">FENE bonds</A>.
+coarse-grained polymer model with <A HREF = "bond_fene.html">FENE bonds</A>.  See
+<A HREF = "#Kremer">(Kremer)</A> for a description of FENE bonds.
 </P>
 <P>The <I>lj/coul</I>, <I>lj</I>, and <I>coul</I> keywords allow the 3 coefficients to
 be set explicitly.  The <I>lj/coul</I> keyword sets both the LJ and
 Coulombic coefficients to the same 3 values.  The <I>lj</I> and <I>coul</I>
 keywords only set either the LJ or Coulombic coefficients.  Use both
 of them if you wish to set the LJ coefficients to different values
 than the Coulombic coefficients.
 </P>
 <P>The <I>angle</I> keyword allows the 1-3 weighting factor to be ignored for
 individual atom pairs if they are not listed as the first and last
 atoms in any angle defined in the simulation or as 1,3 or 2,4 atoms in
 any dihedral defined in the simulation.  For example, imagine the 1-3
 weighting factor is set to 0.5 and you have a linear molecule with 4
 atoms and bonds as follows: 1-2-3-4.  If your data file defines 1-2-3
 as an angle, but does not define 2-3-4 as an angle or 1-2-3-4 as a
 dihedral, then the pairwise interaction between atoms 1 and 3 will
 always be weighted by 0.5, but different force fields use different
 rules for weighting the pairwise interaction between atoms 2 and 4.
 If the <I>angle</I> keyword is specified as <I>yes</I>, then the pairwise
 interaction between atoms 2 and 4 will be unaffected (full weighting
 of 1.0).  If the <I>angle</I> keyword is specified as <I>no</I> which is the
 default, then the 2,4 interaction will also be weighted by 0.5.
 </P>
 <P>The <I>dihedral</I> keyword allows the 1-4 weighting factor to be ignored
 for individual atom pairs if they are not listed as the first and last
 atoms in any dihedral defined in the simulation.  For example, imagine
 the 1-4 weighting factor is set to 0.5 and you have a linear molecule
 with 5 atoms and bonds as follows: 1-2-3-4-5.  If your data file
 defines 1-2-3-4 as a dihedral, but does not define 2-3-4-5 as a
 dihedral, then the pairwise interaction between atoms 1 and 4 will
 always be weighted by 0.5, but different force fields use different
 rules for weighting the pairwise interaction between atoms 2 and 5.
 If the <I>dihedral</I> keyword is specified as <I>yes</I>, then the pairwise
 interaction between atoms 2 and 5 will be unaffected (full weighting
 of 1.0).  If the <I>dihedral</I> keyword is specified as <I>no</I> which is the
 default, then the 2,5 interaction will also be weighted by 0.5.
 </P>
 <P>The <I>extra</I> keyword is used when additional bonds will be created
 during a simulation run, e.g. by the <A HREF = "fix_bond_create.html">fix
 bond/create</A> command.  A list of 1-2,1-3,1-4
 neighbors for each atom is calculated and stored by LAMMPS.  If new
 bonds are created, the list needs to grow.  Using the <I>extra</I> keyword
 leaves empty space in the list for N additional bonds to be added.  If
 you do not do this, you may get an error when bonds are added.
 </P>
 <P><B>Restrictions:</B> none
 </P>
 <P><B>Related commands:</B>
 </P>
 <P><A HREF = "delete_bonds.html">delete_bonds</A>, <A HREF = "fix_bond_create.html">fix bond/create</A>
 </P>
 <P><B>Default:</B>
 </P>
 <P>All 3 Lennard-Jones and 3 Coulobmic weighting coefficients = 0.0,
 angle = no, dihedral = no, and extra = 0.
 </P>
+<HR>
+
+<A NAME = "Cornell"></A>
+
+<P><B>(Cornell)</B> Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
+Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
+</P>
+<A NAME = "Kremer"></A>
+
+<P><B>(Kremer)</B> Kremer, Grest, J Chem Phys, 92, 5057 (1990).
+</P>
+<A NAME = "MacKerell"></A>
+
+<P><B>(MacKerell)</B> MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
+Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
+</P>
+<A NAME = "Mayo"></A>
+
+<P><B>(Mayo)</B> Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
+(1990).
+</P>
 </HTML>
diff --git a/doc/special_bonds.txt b/doc/special_bonds.txt
index e0d2715ae..5988115dc 100644
--- a/doc/special_bonds.txt
+++ b/doc/special_bonds.txt
@@ -1,145 +1,171 @@
 "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
 
 special_bonds command :h3
 
 [Syntax:]
 
 special_bonds keyword values ... :pre
 
 one or more keyword/value pairs may be appended :ulb,l
-keyword = {amber} or {charmm} or {fene} or {lj/coul} or {lj} or {coul} or {angle} or {dihedral} or {extra} :l
+keyword = {amber} or {charmm} or {dreiding} or {fene} or {lj/coul} or {lj} or {coul} or {angle} or {dihedral} or {extra} :l
   {amber} values = none
   {charmm} values = none
+  {dreiding} values = none
   {fene} values = none
   {lj/coul} values = w1,w2,w3
     w1,w2,w3 = weights (0.0 to 1.0) on pairwise Lennard-Jones and Coulombic interactions
   {lj} values = w1,w2,w3
     w1,w2,w3 = weights (0.0 to 1.0) on pairwise Lennard-Jones interactions
   {coul} values = w1,w2,w3
     w1,w2,w3 = weights (0.0 to 1.0) on pairwise Coulombic interactions
   {angle} value = {yes} or {no}
   {dihedral} value = {yes} or {no}
   {extra} value = N
     N = number of extra 1-2,1-3,1-4 interactions to save space for :pre
 :ule
 
 Examples:
 
 special_bonds amber
 special_bonds charmm
 special_bonds fene dihedral no
 special_bonds lj/coul 0.0 0.0 0.5 angle yes dihedral yes
 special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.0 dihedral yes
 special_bonds lj/coul 0 1 1 extra 2 :pre
 
 [Description:]
 
 Set weighting coefficients for pairwise energy and force contributions
 from atom pairs that are also bonded to each other directly or
 indirectly.  For Lennard-Jones (LJ) and Coulombic pairwise
 interactions, these coefficients come in sets of three.  The 1st
 coefficient is the weighting factor on 1-2 atom pairs, which are those
 directly bonded to each other.  The 2nd coefficient is the weighting
 factor on 1-3 atom pairs which are those separated by 2 bonds
 (e.g. the two H atoms in a water molecule).  The 3rd coefficient is
 the weighting factor on 1-4 atom pairs which are those separated by 3
 bonds (e.g. the 1st and 4th atoms in a dihedral interaction).  Thus if
 the 1-2 coefficient is set to 0.0, then the pairwise interaction is
 effectively turned off for all pairs of atoms bonded to each other.
 
 IMPORTANT NOTE: For purposes of computing weighted pairwise
 interactions, 1-3 and 1-4 interactions are not defined from the list
 of angles or dihedrals used by the simulation.  Rather, they are
 inferred topologically from the set of bonds defined when the
 simulation is defined from a data or restart file (see
 "read_data"_read_data.html or "read_restart"_read_restart.html
 commands).  Thus the set of 1-2,1-3,1-4 interactions that the weights
 apply to is the same whether angle and dihedral potentials are
 computed or not, and remains the same even if bonds are constrained,
 or turned off, or removed during a simulation.
 
 The two exceptions to this rule are (a) if the {angle} or {dihedral}
 keywords are set to {yes} (see below), or (b) if the
 "delete_bonds"_delete_bonds.html command is used with the {special}
 option that recomputes the 1-2,1-3,1-4 topologies after bonds are
 deleted; see the "delete_bonds"_delete_bonds.html command for more
 details.
 
 The {amber} keyword sets the 3 coefficients to 0.0, 0.0, 0.5 for LJ
 interactions and to 0.0, 0.0, 0.8333 for Coulombic interactions, which
 is the default for a commonly used version of the AMBER force field,
-where the last value is really 5/6.
+where the last value is really 5/6.  See "(Cornell)"_#Cornell for a
+description of the AMBER force field.
 
 The {charmm} keyword sets the 3 coefficients to 0.0, 0.0, 0.0 for both
 LJ and Coulombic interactions, which is the default for a commonly
 used version of the CHARMM force field.  Note that in pair styles
 {lj/charmm/coul/charmm} and {lj/charmm/coul/long} the 1-4 coefficients
 are defined explicitly, and these pairwise contributions are computed
 as part of the charmm dihedral style - see the
 "pair_coeff"_pair_coeff.html and "dihedral_style"_dihedral_style.html
-commands for more information.
+commands for more information.  See "(MacKerell)"_#MacKerell for a
+description of the CHARMM force field.
+
+The {dreiding} keyword sets the 3 coefficients to 0.0, 0.0, 1.0 for both
+LJ and Coulombic interactions, which is the default for the Dreiding
+force field, as discussed in "(Mayo)"_#Mayo.
 
 The {fene} keyword sets the 3 coefficients to 0.0, 1.0, 1.0 for both
 LJ and Coulombic interactions, which is consistent with a
-coarse-grained polymer model with "FENE bonds"_bond_fene.html.
+coarse-grained polymer model with "FENE bonds"_bond_fene.html.  See
+"(Kremer)"_#Kremer for a description of FENE bonds.
 
 The {lj/coul}, {lj}, and {coul} keywords allow the 3 coefficients to
 be set explicitly.  The {lj/coul} keyword sets both the LJ and
 Coulombic coefficients to the same 3 values.  The {lj} and {coul}
 keywords only set either the LJ or Coulombic coefficients.  Use both
 of them if you wish to set the LJ coefficients to different values
 than the Coulombic coefficients.
 
 The {angle} keyword allows the 1-3 weighting factor to be ignored for
 individual atom pairs if they are not listed as the first and last
 atoms in any angle defined in the simulation or as 1,3 or 2,4 atoms in
 any dihedral defined in the simulation.  For example, imagine the 1-3
 weighting factor is set to 0.5 and you have a linear molecule with 4
 atoms and bonds as follows: 1-2-3-4.  If your data file defines 1-2-3
 as an angle, but does not define 2-3-4 as an angle or 1-2-3-4 as a
 dihedral, then the pairwise interaction between atoms 1 and 3 will
 always be weighted by 0.5, but different force fields use different
 rules for weighting the pairwise interaction between atoms 2 and 4.
 If the {angle} keyword is specified as {yes}, then the pairwise
 interaction between atoms 2 and 4 will be unaffected (full weighting
 of 1.0).  If the {angle} keyword is specified as {no} which is the
 default, then the 2,4 interaction will also be weighted by 0.5.
 
 The {dihedral} keyword allows the 1-4 weighting factor to be ignored
 for individual atom pairs if they are not listed as the first and last
 atoms in any dihedral defined in the simulation.  For example, imagine
 the 1-4 weighting factor is set to 0.5 and you have a linear molecule
 with 5 atoms and bonds as follows: 1-2-3-4-5.  If your data file
 defines 1-2-3-4 as a dihedral, but does not define 2-3-4-5 as a
 dihedral, then the pairwise interaction between atoms 1 and 4 will
 always be weighted by 0.5, but different force fields use different
 rules for weighting the pairwise interaction between atoms 2 and 5.
 If the {dihedral} keyword is specified as {yes}, then the pairwise
 interaction between atoms 2 and 5 will be unaffected (full weighting
 of 1.0).  If the {dihedral} keyword is specified as {no} which is the
 default, then the 2,5 interaction will also be weighted by 0.5.
 
 The {extra} keyword is used when additional bonds will be created
 during a simulation run, e.g. by the "fix
 bond/create"_fix_bond_create.html command.  A list of 1-2,1-3,1-4
 neighbors for each atom is calculated and stored by LAMMPS.  If new
 bonds are created, the list needs to grow.  Using the {extra} keyword
 leaves empty space in the list for N additional bonds to be added.  If
 you do not do this, you may get an error when bonds are added.
 
 [Restrictions:] none
 
 [Related commands:]
 
 "delete_bonds"_delete_bonds.html, "fix bond/create"_fix_bond_create.html
 
 [Default:]
 
 All 3 Lennard-Jones and 3 Coulobmic weighting coefficients = 0.0,
 angle = no, dihedral = no, and extra = 0.
+
+:line
+
+:link(Cornell)
+[(Cornell)] Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
+Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
+
+:link(Kremer)
+[(Kremer)] Kremer, Grest, J Chem Phys, 92, 5057 (1990).
+
+:link(MacKerell)
+[(MacKerell)] MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
+Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
+
+:link(Mayo)
+[(Mayo)] Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
+(1990).
+