git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15293 f3b2605a-c512-4ea7-a41b-209d697bcdaa

doc/html/Manual.html

@@ -135,7 +135,7 @@
   <H1></H1><div class="section" id="lammps-documentation">
 <h1>LAMMPS Documentation</h1>
 <div class="section" id="jul-2016-version">
-<h2>1 Jul 2016 version</h2>
+<h2>9 Jul 2016 version</h2>
 </div>
 <div class="section" id="version-info">
 <h2>Version info:</h2>

doc/html/_sources/Manual.txt

@@ -5,7 +5,7 @@
 LAMMPS Documentation
 ====================
 
-1 Jul 2016 version
+9 Jul 2016 version
 ------------------
 
 Version info:

doc/html/_sources/compute_modify.txt

@@ -19,8 +19,6 @@ Syntax
        N = # of extra degrees of freedom to subtract
      *dynamic* value = *yes* or *no*
        yes/no = do or do not recompute the number of atoms contributing to the temperature
-     *thermo* value = *yes* or *no*
-       yes/no = do or do not add contributions from fixes to the potential energy
 
 
 
@@ -51,18 +49,12 @@ for an example.
 
 The *dynamic* keyword determines whether the number of atoms N in the
 compute group is re-computed each time a temperature is computed.
-Only compute styles that compute a temperature use this option.  By
+Only compute styles that calculate a temperature use this option.  By
 default, N is assumed to be constant.  If you are adding atoms to the
 system (see the :doc:`fix pour <fix_pour>` or :doc:`fix deposit <fix_deposit>` commands) or expect atoms to be lost
 (e.g. due to evaporation), then this option should be used to insure
 the temperature is correctly normalized.
 
-The *thermo* keyword determines whether the potential energy
-contribution calculated by some :doc:`fixes <fix>` is added to the
-potential energy calculated by the compute.  Currently, only the
-compute of style *pe* uses this option.  See the doc pages for
-:doc:`individual fixes <fix>` for details.
-
 Restrictions
 """"""""""""
  none
@@ -76,8 +68,7 @@ Default
 """""""
 
 The option defaults are extra = 2 or 3 for 2d or 3d systems and
-dynamic = no.  Thermo is *yes* if the compute of style *pe* was
+dynamic = no.
-defined with no extra keywords; otherwise it is *no*\ .
 
 
 .. _lws: http://lammps.sandia.gov

doc/html/_sources/compute_pe.txt

@@ -13,7 +13,7 @@ Syntax
 * ID, group-ID are documented in :doc:`compute <compute>` command
 * pe = style name of this compute command
 * zero or more keywords may be appended
-* keyword = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace*
+* keyword = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace* or *fix*
 
 Examples
 """"""""
@@ -35,8 +35,9 @@ via the :doc:`compute reduce <compute_reduce>` command.
 The energy is calculated by the various pair, bond, etc potentials
 defined for the simulation.  If no extra keywords are listed, then the
 potential energy is the sum of pair, bond, angle, dihedral, improper,
-and kspace (long-range) energy.  If any extra keywords are listed,
+kspace (long-range), and fix energy.  I.e. it is as if all the
-then only those components are summed to compute the potential energy.
+keywords were listed.  If any extra keywords are listed, then only
+those components are summed to compute the potential energy.
 
 The Kspace contribution requires 1 extra FFT each timestep the energy
 is calculated, if using the PPPM solver via the :doc:`kspace_style pppm <kspace_style>` command.  Thus it can increase the cost of the
@@ -44,12 +45,14 @@ PPPM calculation if it is needed on a large fraction of the simulation
 timesteps.
 
 Various fixes can contribute to the total potential energy of the
-system.  See the doc pages for :doc:`individual fixes <fix>` for
+system if the *fix* contribution is included.  See the doc pages for
-details.  The *thermo* option of the
+:doc:`individual fixes <fix>` for details of which ones compute a
-:doc:`compute_modify <compute_modify>` command determines whether these
+potential energy.
-contributions are added into the computed potential energy.  If no
+
-keywords are specified the default is *yes*\ .  If any keywords are
+.. note::
-specified, the default is *no*\ .
+
+   The :doc:`fix_modify energy yes <fix_modify>` command must also be
+   specified if a fix is to contribute potential energy to this command.
 
 A compute of this style with the ID of "thermo_pe" is created when
 LAMMPS starts up, as if this command were in the input script:

doc/html/_sources/compute_pe_atom.txt

@@ -13,7 +13,7 @@ Syntax
 * ID, group-ID are documented in :doc:`compute <compute>` command
 * pe/atom = style name of this compute command
 * zero or more keywords may be appended
-* keyword = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace*
+* keyword = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace* or *fix*
 
 Examples
 """"""""
@@ -34,9 +34,9 @@ you want the potential energy of the entire system.
 The per-atom energy is calculated by the various pair, bond, etc
 potentials defined for the simulation.  If no extra keywords are
 listed, then the potential energy is the sum of pair, bond, angle,
-dihedral,improper, and kspace energy.  If any extra keywords are
+dihedral,improper, kspace (long-range), and fix energy.  I.e. it is as
-listed, then only those components are summed to compute the potential
+if all the keywords were listed.  If any extra keywords are listed,
-energy.
+then only those components are summed to compute the potential energy.
 
 Note that the energy of each atom is due to its interaction with all
 other atoms in the simulation, not just with other atoms in the group.
@@ -57,6 +57,17 @@ For PPPM, the calcluation requires 1 extra FFT each timestep that
 per-atom energy is calculated.  Thie `document <PDF/kspace.pdf>`_
 describes how the long-range per-atom energy calculation is performed.
 
+Various fixes can contribute to the per-atom potential energy of the
+system if the *fix* contribution is included.  See the doc pages for
+:doc:`individual fixes <fix>` for details of which ones compute a
+per-atom potential energy.
+
+.. note::
+
+   The :doc:`fix_modify energy yes <fix_modify>` command must also be
+   specified if a fix is to contribute per-atom potential energy to this
+   command.
+
 As an example of per-atom potential energy compared to total potential
 energy, these lines in an input script should yield the same result
 in the last 2 columns of thermo output:

doc/html/_sources/fix_modify.txt

@@ -57,13 +57,21 @@ default method for computing P.
 
 For fixes that calculate a contribution to the potential energy of the
 system, the *energy* keyword will include that contribution in
-thermodynamic output of potential energy.  See the
+thermodynamic output of potential energy.  This is because the *energy
-:doc:`thermo_style <thermo_style>` command for info on how potential
+yes* setting must be specfied to include the fix's global or per-atom
-energy is output.  The contribution by itself can be printed by using
+energy in the calculation performed by the :doc:`compute pe <compute_pe>` or :doc:`compute pe/atom <compute_pe_atom>`
-the keyword f_ID in the thermo_style custom command, where ID is the
+commands.  See the :doc:`thermo_style <thermo_style>` command for info
-fix-ID of the appropriate fix.  Note that you must use this setting
+on how potential energy is output.  For fixes that tally a global
-for a fix if you are using it when performing an :doc:`energy minimization <minimize>` and if you want the energy and forces it
+energy, it can be printed by using the keyword f_ID in the
-produces to be part of the optimization criteria.
+thermo_style custom command, where ID is the fix-ID of the appropriate
+fix.
+
+.. note::
+
+   You must also specify the *energy yes* setting for a fix if you
+   are using it when performing an :doc:`energy minimization <minimize>`
+   and if you want the energy and forces it produces to be part of the
+   optimization criteria.
 
 For fixes that set or modify forces, it may be possible to select at
 which :doc:`r-RESPA <run_style>` level the fix operates via the *respa*

doc/html/_sources/fix_rx.txt

@@ -8,23 +8,30 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group-ID rx file localTemp solver ...
+   fix ID group-ID rx file localTemp matrix solver minSteps ...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * rx = style name of this fix command
 * file = filename containing the reaction kinetic equations and Arrhenius parameters
 * localTemp = *none,lucy* = no local temperature averaging or local temperature defined through Lucy weighting function
-* solver = *lammps_rk4* = Explicit 4th order Runge-Kutta method
+* matrix = *sparse, dense* format for the stoichiometric matrix
-* minSteps = # of steps for rk4 solver
+* solver = *lammps_rk4,rkf45* = rk4 is an explicit 4th order Runge-Kutta method; rkf45 is an adaptive 4th-order Runge-Kutta-Fehlberg method
+* minSteps = # of steps for rk4 solver or minimum # of steps for rkf45 (rk4 or rkf45)
+* maxSteps = maximum number of steps for the rkf45 solver (rkf45 only)
+* relTol = relative tolerance for the rkf45 solver (rkf45 only)
+* absTol = absolute tolernace for the rkf45 solver (rkf45 only)
+* diag   = Diagnostics frequency for the rkf45 solver (optional, rkf45 only)
 
 Examples
 """"""""
 
 .. parsed-literal::
 
-   fix 1 all rx kinetics.rx none lammps_rk4 
+   fix 1 all rx kinetics.rx none dense lammps_rk4 
-   fix 1 all rx kinetics.rx none lammps_rk4 1
+   fix 1 all rx kinetics.rx none sparse lammps_rk4 1
-   fix 1 all rx kinetics.rx lucy lammps_rk4 10
+   fix 1 all rx kinetics.rx lucy sparse lammps_rk4 10 
+   fix 1 all rx kinetics.rx none dense rkf45 1 100 1e-6 1e-8
+   fix 1 all rx kinetics.rx none dense rkf45 1 100 1e-6 1e-8 -1
 
 Description
 """""""""""
@@ -49,13 +56,47 @@ of *m* ordinary differential equations (ODEs) that describe the change
 in concentration of a given species as a function of time are then
 constructed based on the *n* reaction rate equations.
 
-The ODE systems are solved over the full DPD timestep *dt* using a 4th
+The ODE systems are solved over the full DPD timestep *dt* using either a 4th
 order Runge-Kutta *rk4* method with a fixed step-size *h*\ , specified
-by the *lammps_rk4* keyword. The number of ODE steps per DPD timestep
+by the *lammps_rk4* keyword, or a 4th order Runge-Kutta-Fehlberg (rkf45) method 
+with an adaptive step-size for *h*\ . The number of ODE steps per DPD timestep
 for the rk4 method is optionally specified immediately after the rk4
 keyword. The ODE step-size is set as *dt/num_steps*\ . Smaller
 step-sizes tend to yield more accurate results but there is not
-control on the error.
+control on the error. For error control, use the rkf45 ODE solver.
+
+The rkf45 method adjusts the step-size so that the local truncation error is held 
+within the specified absolute and relative tolerances. The initial step-size *h0* 
+can be specified by the user or estimated internally. It is recommeded that the user
+specify *h0* since this will generally reduced the number of ODE integration steps 
+required. *h0* is defined as *dt / min_steps* if min_steps >= 1. If min_steps == 0, 
+*h0* is estimated such that an explicit Euler method would likely produce
+an acceptable solution. This is generally overly conservative for the 4th-order
+method and users are advised to specify *h0* as some fraction of the DPD timestep.
+For small DPD timesteps, only one step may be necessary depending upon the tolerances.
+Note that more than min_steps ODE steps may be taken depending upon the ODE stiffness
+but no more than max_steps will be taken. If max_steps is reached, an error warning
+is printed and the simulation is stopped.
+
+After each ODE step, the solution error *e* is tested and weighted using the absTol
+and relTol values. The error vector is weighted as *e* / (relTol * |\ *u*\ | + absTol) 
+where *u* is the solution vector. If the norm of the error is <= 1, the solution is 
+accepted, *h* is increased by a proportional amount, and the next ODE step is begun.
+Otherwise, *h* is shrunk and the ODE step is repeated.
+
+Run-time diagnostics are available for the rkf45 ODE solver. The frequency 
+(in time-steps) that diagnostics are reported is controlled by the last (optional) 
+12th argument. A negative frequency means that diagnostics are reported once at the 
+end of each run. A positive value N means that the diagnostics are reported once 
+per N time-steps.
+
+The diagnostics report the average # of integrator steps and RHS function evaluations 
+and run-time per ODE as well as the the average/RMS/min/max per process. If the 
+reporting frequency is 1, the RMS/min/max per ODE are also reported. The per ODE 
+statistics can be used to adjust the tolerance and min/max step parameters. The 
+statistics per MPI process can be useful to examine any load imbalance caused by the 
+adaptive ODE solver. (Some DPD particles can take longer to solve than others. This 
+can lead to an imbalance across the MPI processes.)
 
 
 ----------
@@ -100,6 +141,17 @@ where the Lucy function is expressed as:
 
 The self-particle interaction is included in the above equation.
 
+The stoichiometric coefficients for the reaction mechanism are stored
+in either a sparse or dense matrix format. The dense matrix should only be
+used for small reaction mechanisms. The sparse matrix should be used when there
+are many reactions (e.g., more than 5). This allows the number of reactions and
+species to grow while keeping the computational cost tractable. The matrix
+format can be specified as using either the *sparse* or *dense* keywords.
+If all stoichiometric coefficients for a reaction are small integers (whole
+numbers <= 3), a fast exponential function is used. This can save significant
+computational time so users are encouraged to use integer coefficients
+where possible.
+
 
 ----------
 

doc/html/compute_modify.html

@@ -141,8 +141,6 @@
   N = # of extra degrees of freedom to subtract
 <em>dynamic</em> value = <em>yes</em> or <em>no</em>
   yes/no = do or do not recompute the number of atoms contributing to the temperature
-<em>thermo</em> value = <em>yes</em> or <em>no</em>
-  yes/no = do or do not add contributions from fixes to the potential energy
 </pre>
 </div>
 <div class="section" id="examples">
@@ -169,16 +167,11 @@ the <em>extra</em> parameter if you need to add degrees-of-freedom.</p>
 <p>for an example.</p>
 <p>The <em>dynamic</em> keyword determines whether the number of atoms N in the
 compute group is re-computed each time a temperature is computed.
-Only compute styles that compute a temperature use this option.  By
+Only compute styles that calculate a temperature use this option.  By
 default, N is assumed to be constant.  If you are adding atoms to the
 system (see the <a class="reference internal" href="fix_pour.html"><span class="doc">fix pour</span></a> or <a class="reference internal" href="fix_deposit.html"><span class="doc">fix deposit</span></a> commands) or expect atoms to be lost
 (e.g. due to evaporation), then this option should be used to insure
 the temperature is correctly normalized.</p>
-<p>The <em>thermo</em> keyword determines whether the potential energy
-contribution calculated by some <a class="reference internal" href="fix.html"><span class="doc">fixes</span></a> is added to the
-potential energy calculated by the compute.  Currently, only the
-compute of style <em>pe</em> uses this option.  See the doc pages for
-<a class="reference internal" href="fix.html"><span class="doc">individual fixes</span></a> for details.</p>
 </div>
 <div class="section" id="restrictions">
 <h2>Restrictions</h2>
@@ -192,8 +185,7 @@ compute of style <em>pe</em> uses this option.  See the doc pages for
 <div class="section" id="default">
 <h2>Default</h2>
 <p>The option defaults are extra = 2 or 3 for 2d or 3d systems and
-dynamic = no.  Thermo is <em>yes</em> if the compute of style <em>pe</em> was
+dynamic = no.</p>
-defined with no extra keywords; otherwise it is <em>no</em>.</p>
 </div>
 </div>
 

doc/html/compute_pe.html

@@ -135,7 +135,7 @@
 <li>ID, group-ID are documented in <a class="reference internal" href="compute.html"><span class="doc">compute</span></a> command</li>
 <li>pe = style name of this compute command</li>
 <li>zero or more keywords may be appended</li>
-<li>keyword = <em>pair</em> or <em>bond</em> or <em>angle</em> or <em>dihedral</em> or <em>improper</em> or <em>kspace</em></li>
+<li>keyword = <em>pair</em> or <em>bond</em> or <em>angle</em> or <em>dihedral</em> or <em>improper</em> or <em>kspace</em> or <em>fix</em></li>
 </ul>
 </div>
 <div class="section" id="examples">
@@ -155,19 +155,22 @@ via the <a class="reference internal" href="compute_reduce.html"><span class="do
 <p>The energy is calculated by the various pair, bond, etc potentials
 defined for the simulation.  If no extra keywords are listed, then the
 potential energy is the sum of pair, bond, angle, dihedral, improper,
-and kspace (long-range) energy.  If any extra keywords are listed,
+kspace (long-range), and fix energy.  I.e. it is as if all the
-then only those components are summed to compute the potential energy.</p>
+keywords were listed.  If any extra keywords are listed, then only
+those components are summed to compute the potential energy.</p>
 <p>The Kspace contribution requires 1 extra FFT each timestep the energy
 is calculated, if using the PPPM solver via the <a class="reference internal" href="kspace_style.html"><span class="doc">kspace_style pppm</span></a> command.  Thus it can increase the cost of the
 PPPM calculation if it is needed on a large fraction of the simulation
 timesteps.</p>
 <p>Various fixes can contribute to the total potential energy of the
-system.  See the doc pages for <a class="reference internal" href="fix.html"><span class="doc">individual fixes</span></a> for
+system if the <em>fix</em> contribution is included.  See the doc pages for
-details.  The <em>thermo</em> option of the
+<a class="reference internal" href="fix.html"><span class="doc">individual fixes</span></a> for details of which ones compute a
-<a class="reference internal" href="compute_modify.html"><span class="doc">compute_modify</span></a> command determines whether these
+potential energy.</p>
-contributions are added into the computed potential energy.  If no
+<div class="admonition note">
-keywords are specified the default is <em>yes</em>.  If any keywords are
+<p class="first admonition-title">Note</p>
-specified, the default is <em>no</em>.</p>
+<p class="last">The <a class="reference internal" href="fix_modify.html"><span class="doc">fix_modify energy yes</span></a> command must also be
+specified if a fix is to contribute potential energy to this command.</p>
+</div>
 <p>A compute of this style with the ID of &#8220;thermo_pe&#8221; is created when
 LAMMPS starts up, as if this command were in the input script:</p>
 <div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">thermo_pe</span> <span class="nb">all</span> <span class="n">pe</span>

doc/html/compute_pe_atom.html

@@ -135,7 +135,7 @@
 <li>ID, group-ID are documented in <a class="reference internal" href="compute.html"><span class="doc">compute</span></a> command</li>
 <li>pe/atom = style name of this compute command</li>
 <li>zero or more keywords may be appended</li>
-<li>keyword = <em>pair</em> or <em>bond</em> or <em>angle</em> or <em>dihedral</em> or <em>improper</em> or <em>kspace</em></li>
+<li>keyword = <em>pair</em> or <em>bond</em> or <em>angle</em> or <em>dihedral</em> or <em>improper</em> or <em>kspace</em> or <em>fix</em></li>
 </ul>
 </div>
 <div class="section" id="examples">
@@ -154,9 +154,9 @@ you want the potential energy of the entire system.</p>
 <p>The per-atom energy is calculated by the various pair, bond, etc
 potentials defined for the simulation.  If no extra keywords are
 listed, then the potential energy is the sum of pair, bond, angle,
-dihedral,improper, and kspace energy.  If any extra keywords are
+dihedral,improper, kspace (long-range), and fix energy.  I.e. it is as
-listed, then only those components are summed to compute the potential
+if all the keywords were listed.  If any extra keywords are listed,
-energy.</p>
+then only those components are summed to compute the potential energy.</p>
 <p>Note that the energy of each atom is due to its interaction with all
 other atoms in the simulation, not just with other atoms in the group.</p>
 <p>For an energy contribution produced by a small set of atoms (e.g. 4
@@ -172,6 +172,16 @@ as specified by the <a class="reference internal" href="kspace_style.html"><span
 For PPPM, the calcluation requires 1 extra FFT each timestep that
 per-atom energy is calculated.  Thie <a class="reference external" href="PDF/kspace.pdf">document</a>
 describes how the long-range per-atom energy calculation is performed.</p>
+<p>Various fixes can contribute to the per-atom potential energy of the
+system if the <em>fix</em> contribution is included.  See the doc pages for
+<a class="reference internal" href="fix.html"><span class="doc">individual fixes</span></a> for details of which ones compute a
+per-atom potential energy.</p>
+<div class="admonition note">
+<p class="first admonition-title">Note</p>
+<p class="last">The <a class="reference internal" href="fix_modify.html"><span class="doc">fix_modify energy yes</span></a> command must also be
+specified if a fix is to contribute per-atom potential energy to this
+command.</p>
+</div>
 <p>As an example of per-atom potential energy compared to total potential
 energy, these lines in an input script should yield the same result
 in the last 2 columns of thermo output:</p>

doc/html/fix_modify.html

@@ -173,13 +173,21 @@ documentation.  Thus this option allows the user to override the
 default method for computing P.</p>
 <p>For fixes that calculate a contribution to the potential energy of the
 system, the <em>energy</em> keyword will include that contribution in
-thermodynamic output of potential energy.  See the
+thermodynamic output of potential energy.  This is because the <em>energy
-<a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style</span></a> command for info on how potential
+yes</em> setting must be specfied to include the fix&#8217;s global or per-atom
-energy is output.  The contribution by itself can be printed by using
+energy in the calculation performed by the <a class="reference internal" href="compute_pe.html"><span class="doc">compute pe</span></a> or <a class="reference internal" href="compute_pe_atom.html"><span class="doc">compute pe/atom</span></a>
-the keyword f_ID in the thermo_style custom command, where ID is the
+commands.  See the <a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style</span></a> command for info
-fix-ID of the appropriate fix.  Note that you must use this setting
+on how potential energy is output.  For fixes that tally a global
-for a fix if you are using it when performing an <a class="reference internal" href="minimize.html"><span class="doc">energy minimization</span></a> and if you want the energy and forces it
+energy, it can be printed by using the keyword f_ID in the
-produces to be part of the optimization criteria.</p>
+thermo_style custom command, where ID is the fix-ID of the appropriate
+fix.</p>
+<div class="admonition note">
+<p class="first admonition-title">Note</p>
+<p class="last">You must also specify the <em>energy yes</em> setting for a fix if you
+are using it when performing an <a class="reference internal" href="minimize.html"><span class="doc">energy minimization</span></a>
+and if you want the energy and forces it produces to be part of the
+optimization criteria.</p>
+</div>
 <p>For fixes that set or modify forces, it may be possible to select at
 which <a class="reference internal" href="run_style.html"><span class="doc">r-RESPA</span></a> level the fix operates via the <em>respa</em>
 keyword. The RESPA level at which the fix is active can be selected.

doc/html/fix_rx.html

@@ -128,7 +128,7 @@
 <span id="index-0"></span><h1>fix rx command</h1>
 <div class="section" id="syntax">
 <h2>Syntax</h2>
-<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="n">ID</span> <span class="n">group</span><span class="o">-</span><span class="n">ID</span> <span class="n">rx</span> <span class="n">file</span> <span class="n">localTemp</span> <span class="n">solver</span> <span class="o">...</span>
+<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="n">ID</span> <span class="n">group</span><span class="o">-</span><span class="n">ID</span> <span class="n">rx</span> <span class="n">file</span> <span class="n">localTemp</span> <span class="n">matrix</span> <span class="n">solver</span> <span class="n">minSteps</span> <span class="o">...</span>
 </pre></div>
 </div>
 <ul class="simple">
@@ -136,15 +136,22 @@
 <li>rx = style name of this fix command</li>
 <li>file = filename containing the reaction kinetic equations and Arrhenius parameters</li>
 <li>localTemp = <em>none,lucy</em> = no local temperature averaging or local temperature defined through Lucy weighting function</li>
-<li>solver = <em>lammps_rk4</em> = Explicit 4th order Runge-Kutta method</li>
+<li>matrix = <em>sparse, dense</em> format for the stoichiometric matrix</li>
-<li>minSteps = # of steps for rk4 solver</li>
+<li>solver = <em>lammps_rk4,rkf45</em> = rk4 is an explicit 4th order Runge-Kutta method; rkf45 is an adaptive 4th-order Runge-Kutta-Fehlberg method</li>
+<li>minSteps = # of steps for rk4 solver or minimum # of steps for rkf45 (rk4 or rkf45)</li>
+<li>maxSteps = maximum number of steps for the rkf45 solver (rkf45 only)</li>
+<li>relTol = relative tolerance for the rkf45 solver (rkf45 only)</li>
+<li>absTol = absolute tolernace for the rkf45 solver (rkf45 only)</li>
+<li>diag   = Diagnostics frequency for the rkf45 solver (optional, rkf45 only)</li>
 </ul>
 </div>
 <div class="section" id="examples">
 <h2>Examples</h2>
-<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">none</span> <span class="n">lammps_rk4</span>
+<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">none</span> <span class="n">dense</span> <span class="n">lammps_rk4</span>
-<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">none</span> <span class="n">lammps_rk4</span> <span class="mi">1</span>
+<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">none</span> <span class="n">sparse</span> <span class="n">lammps_rk4</span> <span class="mi">1</span>
-<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">lucy</span> <span class="n">lammps_rk4</span> <span class="mi">10</span>
+<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">lucy</span> <span class="n">sparse</span> <span class="n">lammps_rk4</span> <span class="mi">10</span>
+<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">none</span> <span class="n">dense</span> <span class="n">rkf45</span> <span class="mi">1</span> <span class="mi">100</span> <span class="mi">1</span><span class="n">e</span><span class="o">-</span><span class="mi">6</span> <span class="mi">1</span><span class="n">e</span><span class="o">-</span><span class="mi">8</span>
+<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">rx</span> <span class="n">kinetics</span><span class="o">.</span><span class="n">rx</span> <span class="n">none</span> <span class="n">dense</span> <span class="n">rkf45</span> <span class="mi">1</span> <span class="mi">100</span> <span class="mi">1</span><span class="n">e</span><span class="o">-</span><span class="mi">6</span> <span class="mi">1</span><span class="n">e</span><span class="o">-</span><span class="mi">8</span> <span class="o">-</span><span class="mi">1</span>
 </pre></div>
 </div>
 </div>
@@ -162,13 +169,43 @@ of <em>n</em> reaction equations will contain a total of <em>m</em> species.  A
 of <em>m</em> ordinary differential equations (ODEs) that describe the change
 in concentration of a given species as a function of time are then
 constructed based on the <em>n</em> reaction rate equations.</p>
-<p>The ODE systems are solved over the full DPD timestep <em>dt</em> using a 4th
+<p>The ODE systems are solved over the full DPD timestep <em>dt</em> using either a 4th
 order Runge-Kutta <em>rk4</em> method with a fixed step-size <em>h</em>, specified
-by the <em>lammps_rk4</em> keyword. The number of ODE steps per DPD timestep
+by the <em>lammps_rk4</em> keyword, or a 4th order Runge-Kutta-Fehlberg (rkf45) method
+with an adaptive step-size for <em>h</em>. The number of ODE steps per DPD timestep
 for the rk4 method is optionally specified immediately after the rk4
 keyword. The ODE step-size is set as <em>dt/num_steps</em>. Smaller
 step-sizes tend to yield more accurate results but there is not
-control on the error.</p>
+control on the error. For error control, use the rkf45 ODE solver.</p>
+<p>The rkf45 method adjusts the step-size so that the local truncation error is held
+within the specified absolute and relative tolerances. The initial step-size <em>h0</em>
+can be specified by the user or estimated internally. It is recommeded that the user
+specify <em>h0</em> since this will generally reduced the number of ODE integration steps
+required. <em>h0</em> is defined as <em>dt / min_steps</em> if min_steps &gt;= 1. If min_steps == 0,
+<em>h0</em> is estimated such that an explicit Euler method would likely produce
+an acceptable solution. This is generally overly conservative for the 4th-order
+method and users are advised to specify <em>h0</em> as some fraction of the DPD timestep.
+For small DPD timesteps, only one step may be necessary depending upon the tolerances.
+Note that more than min_steps ODE steps may be taken depending upon the ODE stiffness
+but no more than max_steps will be taken. If max_steps is reached, an error warning
+is printed and the simulation is stopped.</p>
+<p>After each ODE step, the solution error <em>e</em> is tested and weighted using the absTol
+and relTol values. The error vector is weighted as <em>e</em> / (relTol * <a href="#id1"><span class="problematic" id="id2">|</span></a><em>u</em>| + absTol)
+where <em>u</em> is the solution vector. If the norm of the error is &lt;= 1, the solution is
+accepted, <em>h</em> is increased by a proportional amount, and the next ODE step is begun.
+Otherwise, <em>h</em> is shrunk and the ODE step is repeated.</p>
+<p>Run-time diagnostics are available for the rkf45 ODE solver. The frequency
+(in time-steps) that diagnostics are reported is controlled by the last (optional)
+12th argument. A negative frequency means that diagnostics are reported once at the
+end of each run. A positive value N means that the diagnostics are reported once
+per N time-steps.</p>
+<p>The diagnostics report the average # of integrator steps and RHS function evaluations
+and run-time per ODE as well as the the average/RMS/min/max per process. If the
+reporting frequency is 1, the RMS/min/max per ODE are also reported. The per ODE
+statistics can be used to adjust the tolerance and min/max step parameters. The
+statistics per MPI process can be useful to examine any load imbalance caused by the
+adaptive ODE solver. (Some DPD particles can take longer to solve than others. This
+can lead to an imbalance across the MPI processes.)</p>
 <hr class="docutils" />
 <p>The filename specifies a file that contains the entire set of reaction
 kinetic equations and corresponding Arrhenius parameters.  The format of
@@ -197,6 +234,16 @@ particle internal temperature is defined as:</p>
 <p>where the Lucy function is expressed as:</p>
 <img alt="_images/fix_rx_localTemp2.jpg" class="align-center" src="_images/fix_rx_localTemp2.jpg" />
 <p>The self-particle interaction is included in the above equation.</p>
+<p>The stoichiometric coefficients for the reaction mechanism are stored
+in either a sparse or dense matrix format. The dense matrix should only be
+used for small reaction mechanisms. The sparse matrix should be used when there
+are many reactions (e.g., more than 5). This allows the number of reactions and
+species to grow while keeping the computational cost tractable. The matrix
+format can be specified as using either the <em>sparse</em> or <em>dense</em> keywords.
+If all stoichiometric coefficients for a reaction are small integers (whole
+numbers &lt;= 3), a fast exponential function is used. This can save significant
+computational time so users are encouraged to use integer coefficients
+where possible.</p>
 <hr class="docutils" />
 <p>The format of a tabulated file is as follows (without the
 parenthesized comments):</p>
@@ -224,7 +271,7 @@ temperature exponent <em>n</em>, and the activation energy <em>Ea</em>.</p>
 are used by the <a class="reference internal" href="fix_eos_table_rx.html"><span class="doc">fix eos/table/rx</span></a> command to
 define the thermodynamic properties of each species.  Furthermore, the
 number of species molecules (i.e., concentration) can be specified
-either with the <a class="reference internal" href="set.html"><span class="doc">set</span></a> command using the &#8220;<a href="#id1"><span class="problematic" id="id2">d_</span></a>&#8221; prefix or by
+either with the <a class="reference internal" href="set.html"><span class="doc">set</span></a> command using the &#8220;<a href="#id3"><span class="problematic" id="id4">d_</span></a>&#8221; prefix or by
 reading directly the concentrations from a data file.  For the latter
 case, the <a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> command with the fix keyword
 should be specified, where the fix-ID will be the &#8220;fix rx`ID with a &lt;SPECIES&#8221;&gt;`_ suffix, e.g.</p>