convert pair styles gauss to gromacs

doc/src/Eqs/pair_gauss.tex

@@ -1,9 +0,0 @@
-\documentclass[12pt]{article}
-\pagestyle{empty}
-\begin{document}
-
-$$
-   E = - A \exp(-B r^2) \qquad r < r_c
-$$
-
-\end{document}

doc/src/Eqs/pair_gauss_cut.tex

@@ -1,8 +0,0 @@
-\documentclass[12pt]{article}
-\pagestyle{empty}
-\begin{document}
-
-\begin{eqnarray*}
- E = & \frac{H}{\sigma_h\sqrt{2\pi}} \exp\left[-\frac{(r-r_{mh})^2}{2\sigma_h^2}\right]
-\end{eqnarray*}
-\end{document}

doc/src/Eqs/pair_gayberne.tex

@@ -1,14 +0,0 @@
-\documentstyle[12pt]{article}
-
-\begin{document}
-
-$$ U ( \mathbf{A}_1, \mathbf{A}_2, \mathbf{r}_{12} ) = U_r (
-\mathbf{A}_1, \mathbf{A}_2, \mathbf{r}_{12}, \gamma ) \cdot \eta_{12} (
-\mathbf{A}_1, \mathbf{A}_2, \upsilon ) \cdot \chi_{12} ( \mathbf{A}_1,
-\mathbf{A}_2, \mathbf{r}_{12}, \mu ) $$
-
-$$ U_r = 4 \epsilon ( \varrho^{12} - \varrho^6) $$
-
-$$ \varrho = \frac{\sigma}{ h_{12} + \gamma \sigma} $$
-
-\end{document}

doc/src/Eqs/pair_gayberne2.tex

@@ -1,9 +0,0 @@
-\documentstyle[12pt]{article}
-
-\begin{document}
-
-$$ \epsilon_a = \sigma \cdot { \frac{a}{ b \cdot c } }; \epsilon_b =
-\sigma \cdot { \frac{b}{ a \cdot c } }; \epsilon_c = \sigma \cdot {
-\frac{c}{ a \cdot b } } $$
-
-\end{document}

doc/src/Eqs/pair_gran_hertz.tex

@@ -1,14 +0,0 @@
-\documentclass[12pt]{article}
-
-\begin{document}
-
-$$
-   F_{hz} = \sqrt{\delta} \sqrt{\frac{R_i R_j}{R_i + R_j}} F_{hk} = 
-     \sqrt{\delta} \sqrt{\frac{R_i R_j}{R_i + R_j}} 
-     \Big[ (k_n \delta \mathbf{n}_{ij} -  
-       m_{\mbox{\scriptsize{eff}}} \: \gamma_n \mathbf{ v}_n) -
-       (k_t \mathbf{ \Delta s}_t +
-       m_{\mbox{\scriptsize{eff}}} \: \gamma_t \mathbf{v}_t) \Big]
-$$
-
-\end{document}

doc/src/Eqs/pair_gran_hooke.tex

@@ -1,12 +0,0 @@
-\documentclass[12pt]{article}
-
-\begin{document}
-
-$$
-   F_{hk} = (k_n \delta \mathbf{n}_{ij} -  
-   m_{\mbox{\scriptsize{eff}}} \gamma_n\mathbf{ v}_n) - 
-   (k_t \mathbf{ \Delta s}_t +
-   m_{\mbox{\scriptsize{eff}}} \gamma_t \mathbf{v}_t)
-$$
-
-\end{document}

doc/src/Eqs/pair_gromacs.tex

@@ -1,17 +0,0 @@
-\documentstyle[12pt]{article}
-
-\begin{document}
-
-\begin{eqnarray*}
-E_{LJ} & = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - 
-                       \left(\frac{\sigma}{r}\right)^6 \right] + S_{LJ}(r)
-                       \qquad r < r_c \\
-E_C & = & \frac{C q_i q_j}{\epsilon  r} + S_C(r) \qquad r < r_c \\
-S(r) & = & C \qquad r < r_1 \\
-S(r) & = & \frac{A}{3} (r - r_1)^3 + \frac{B}{4} (r - r_1)^4 + C \qquad  r_1 < r < r_c \\
-A & = & (-3 E'(r_c) + (r_c - r_1) E''(r_c))/(r_c - r_1)^2 \\
-B & = & (2 E'(r_c) - (r_c - r_1) E''(r_c))/(r_c - r_1)^3 \\
-C & = & -E(r_c) + \frac{1}{2} (r_c - r_1) E'(r_c) - \frac{1}{12} (r_c - r_1)^2 E''(r_c)) \\
-\end{eqnarray*}                           
-
-\end{document}

doc/src/pair_gauss.rst

@@ -44,11 +44,13 @@ Description
 
 Style *gauss* computes a tethering potential of the form
 
-.. image:: Eqs/pair_gauss.jpg
-   :align: center
+.. math::
+
+   E = - A \exp(-B r^2) \qquad r < r_c
+
 
 between an atom and its corresponding tether site which will typically
-be a frozen atom in the simulation.  Rc is the cutoff.
+be a frozen atom in the simulation.  :math:`r_c` is the cutoff.
 
 The following coefficients must be defined for each pair of atom types
 via the :doc:`pair_coeff <pair_coeff>` command as in the examples above,
@@ -66,17 +68,20 @@ is used.
 Style *gauss/cut* computes a generalized Gaussian interaction potential
 between pairs of particles:
 
-.. image:: Eqs/pair_gauss_cut.jpg
-   :align: center
+.. math::
 
-where H determines together with the standard deviation sigma\_h the
-peak height of the Gaussian function, and r\_mh the peak position.
-Examples of the use of the Gaussian potentials include implicit
-solvent simulations of salt ions :ref:`(Lenart) <Lenart2>` and of surfactants
-:ref:`(Jusufi) <Jusufi2>`.  In these instances the Gaussian potential mimics
-the hydration barrier between a pair of particles. The hydration
-barrier is located at r\_mh and has a width of sigma\_h. The prefactor
-determines the height of the potential barrier.
+   E = \frac{H}{\sigma_h\sqrt{2\pi}} \exp\left[-\frac{(r-r_{mh})^2}{2\sigma_h^2}\right]
+
+
+where H determines together with the standard deviation :math:`\sigma_h`
+the peak height of the Gaussian function, and :math:`r_{mh}` the peak
+position.  Examples of the use of the Gaussian potentials include
+implicit solvent simulations of salt ions :ref:`(Lenart) <Lenart2>` and
+of surfactants :ref:`(Jusufi) <Jusufi2>`.  In these instances the
+Gaussian potential mimics the hydration barrier between a pair of
+particles. The hydration barrier is located at :math:`r_{mh}` and has a
+width of :math:`\sigma_h`. The prefactor determines the height of the
+potential barrier.
 
 The following coefficients must be defined for each pair of atom types
 via the :doc:`pair_coeff <pair_coeff>` command as in the example above,
@@ -85,17 +90,15 @@ or in the data file or restart files read by the
 commands:
 
 * H (energy \* distance units)
-* r\_mh (distance units)
-* sigma\_h (distance units)
+* :math:`r_{mh}` (distance units)
+* :math:`\sigma_h` (distance units)
 * cutoff (distance units)
 
 The last coefficient is optional. If not specified, the global cutoff
 is used.
 
-
 ----------
 
-
 Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are
 functionally the same as the corresponding style without the suffix.
 They have been optimized to run faster, depending on your available

doc/src/pair_gayberne.rst

@@ -41,24 +41,33 @@ The *gayberne* styles compute a Gay-Berne anisotropic LJ interaction
 :ref:`(Berardi) <Berardi>` between pairs of ellipsoidal particles or an
 ellipsoidal and spherical particle via the formulas
 
-.. image:: Eqs/pair_gayberne.jpg
-   :align: center
+.. math::
 
-where A1 and A2 are the transformation matrices from the simulation
-box frame to the body frame and r12 is the center to center vector
-between the particles.  Ur controls the shifted distance dependent
-interaction based on the distance of closest approach of the two
-particles (h12) and the user-specified shift parameter gamma.  When
-both particles are spherical, the formula reduces to the usual
-Lennard-Jones interaction (see details below for when Gay-Berne treats
-a particle as "spherical").
+   U ( \mathbf{A}_1, \mathbf{A}_2, \mathbf{r}_{12} ) = & U_r (
+   \mathbf{A}_1, \mathbf{A}_2, \mathbf{r}_{12}, \gamma ) \cdot \eta_{12} (
+   \mathbf{A}_1, \mathbf{A}_2, \upsilon ) \cdot \chi_{12} ( \mathbf{A}_1,
+   \mathbf{A}_2, \mathbf{r}_{12}, \mu ) \\
+   U_r = & 4 \epsilon ( \varrho^{12} - \varrho^6) \\
+   \varrho = & \frac{\sigma}{ h_{12} + \gamma \sigma}
+
+where A1 and A2 are the transformation matrices from the simulation box
+frame to the body frame and :math:`r_{12}` is the center to center
+vector between the particles.  :math:`U_r` controls the shifted distance
+dependent interaction based on the distance of closest approach of the
+two particles (:math:`h_{12}`) and the user-specified shift parameter
+gamma.  When both particles are spherical, the formula reduces to the
+usual Lennard-Jones interaction (see details below for when Gay-Berne
+treats a particle as "spherical").
 
 For large uniform molecules it has been shown that the energy
 parameters are approximately representable in terms of local contact
 curvatures :ref:`(Everaers) <Everaers2>`:
 
-.. image:: Eqs/pair_gayberne2.jpg
-   :align: center
+.. math::
+
+   \epsilon_a = \sigma \cdot { \frac{a}{ b \cdot c } }; \epsilon_b =
+   \sigma \cdot { \frac{b}{ a \cdot c } }; \epsilon_c = \sigma \cdot {
+   \frac{c}{ a \cdot b } }
 
 The variable names utilized as potential parameters are for the most
 part taken from :ref:`(Everaers) <Everaers2>` in order to be consistent with
@@ -80,69 +89,72 @@ above, or in the data file or restart files read by the
 :doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
 commands, or by mixing as described below:
 
-* epsilon = well depth (energy units)
-* sigma = minimum effective particle radii (distance units)
-* epsilon\_i_a = relative well depth of type I for side-to-side interactions
-* epsilon\_i_b = relative well depth of type I for face-to-face interactions
-* epsilon\_i_c = relative well depth of type I for end-to-end interactions
-* epsilon\_j_a = relative well depth of type J for side-to-side interactions
-* epsilon\_j_b = relative well depth of type J for face-to-face interactions
-* epsilon\_j_c = relative well depth of type J for end-to-end interactions
+* :math:`\epsilon` = well depth (energy units)
+* :math:`\sigma` = minimum effective particle radii (distance units)
+* :math:`\epsilon_{i,a}` = relative well depth of type I for side-to-side interactions
+* :math:`\epsilon_{i,b}` = relative well depth of type I for face-to-face interactions
+* :math:`\epsilon_{i,c}` = relative well depth of type I for end-to-end interactions
+* :math:`\epsilon_{j,a}` = relative well depth of type J for side-to-side interactions
+* :math:`\epsilon_{j,b}` = relative well depth of type J for face-to-face interactions
+* :math:`\epsilon_{j,c}` = relative well depth of type J for end-to-end interactions
 * cutoff (distance units)
 
 The last coefficient is optional.  If not specified, the global
 cutoff specified in the pair\_style command is used.
 
-It is typical with the Gay-Berne potential to define *sigma* as the
-minimum of the 3 shape diameters of the particles involved in an I,I
-interaction, though this is not required.  Note that this is a
-different meaning for *sigma* than the :doc:`pair_style resquared <pair_resquared>` potential uses.
+It is typical with the Gay-Berne potential to define :math:`\sigma` as
+the minimum of the 3 shape diameters of the particles involved in an I,I
+interaction, though this is not required.  Note that this is a different
+meaning for :math:`\sigma` than the :doc:`pair_style resquared
+<pair_resquared>` potential uses.
 
-The epsilon\_i and epsilon\_j coefficients are actually defined for atom
-types, not for pairs of atom types.  Thus, in a series of pair\_coeff
-commands, they only need to be specified once for each atom type.
+The :math:`\epsilon_i` and :math:`\epsilon_j` coefficients are actually
+defined for atom types, not for pairs of atom types.  Thus, in a series
+of pair\_coeff commands, they only need to be specified once for each
+atom type.
 
-Specifically, if any of epsilon\_i_a, epsilon\_i_b, epsilon\_i_c are
-non-zero, the three values are assigned to atom type I.  If all the
-epsilon\_i values are zero, they are ignored.  If any of epsilon\_j_a,
-epsilon\_j_b, epsilon\_j_c are non-zero, the three values are assigned
-to atom type J.  If all three epsilon\_j values are zero, they are
-ignored.  Thus the typical way to define the epsilon\_i and epsilon\_j
-coefficients is to list their values in "pair\_coeff I J" commands when
-I = J, but set them to 0.0 when I != J.  If you do list them when I !=
-J, you should insure they are consistent with their values in other
-pair\_coeff commands, since only the last setting will be in effect.
+Specifically, if any of :math:`\epsilon_{i,a}`, :math:`\epsilon_{i,b}`,
+:math:`\epsilon_{i,c}` are non-zero, the three values are assigned to
+atom type I.  If all the :math:`\epsilon_i` values are zero, they are
+ignored.  If any of :math:`\epsilon_{j,a}`, :math:`\epsilon_{j,b}`,
+:math:`\epsilon_{j,c}` are non-zero, the three values are assigned to
+atom type J.  If all three epsilon\_j values are zero, they are ignored.
+Thus the typical way to define the :math:`\epsilon_i` and
+:math:`\epsilon_j` coefficients is to list their values in "pair\_coeff
+I J" commands when I = J, but set them to 0.0 when I != J.  If you do
+list them when I != J, you should insure they are consistent with their
+values in other pair\_coeff commands, since only the last setting will
+be in effect.
 
 Note that if this potential is being used as a sub-style of
 :doc:`pair_style hybrid <pair_hybrid>`, and there is no "pair\_coeff I I"
 setting made for Gay-Berne for a particular type I (because I-I
 interactions are computed by another hybrid pair potential), then you
-still need to insure the epsilon a,b,c coefficients are assigned to
+still need to insure the :math:`\epsilon` a,b,c coefficients are assigned to
 that type. e.g. in a "pair\_coeff I J" command.
 
 .. note::
 
-   If the epsilon a = b = c for an atom type, and if the shape of
-   the particle itself is spherical, meaning its 3 shape parameters are
-   all the same, then the particle is treated as an LJ sphere by the
+   If the :math:`\epsilon` a = b = c for an atom type, and if the shape
+   of the particle itself is spherical, meaning its 3 shape parameters
+   are all the same, then the particle is treated as an LJ sphere by the
    Gay-Berne potential.  This is significant because if two LJ spheres
    interact, then the simple Lennard-Jones formula is used to compute
    their interaction energy/force using the specified epsilon and sigma
    as the standard LJ parameters.  This is much cheaper to compute than
-   the full Gay-Berne formula.  To treat the particle as a LJ sphere with
-   sigma = D, you should normally set epsilon a = b = c = 1.0, set the
-   pair\_coeff sigma = D, and also set the 3 shape parameters for the
-   particle to D.  The one exception is that if the 3 shape parameters
-   are set to 0.0, which is a valid way in LAMMPS to specify a point
-   particle, then the Gay-Berne potential will treat that as shape
-   parameters of 1.0 (i.e. a LJ particle with sigma = 1), since it
-   requires finite-size particles.  In this case you should still set the
-   pair\_coeff sigma to 1.0 as well.
-
+   the full Gay-Berne formula.  To treat the particle as a LJ sphere
+   with sigma = D, you should normally set :math:`\epsilon` a = b = c =
+   1.0, set the pair\_coeff :math:`\sigma = D`, and also set the 3 shape
+   parameters for the particle to D.  The one exception is that if the 3
+   shape parameters are set to 0.0, which is a valid way in LAMMPS to
+   specify a point particle, then the Gay-Berne potential will treat
+   that as shape parameters of 1.0 (i.e. a LJ particle with
+   :math:`\sigma = 1`), since it requires finite-size particles.  In
+   this case you should still set the pair\_coeff :math:`\sigma` to 1.0
+   as well.
 
 ----------
 
-
 Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are
 functionally the same as the corresponding style without the suffix.
 They have been optimized to run faster, depending on your available

doc/src/pair_gran.rst

@@ -70,13 +70,25 @@ no force between the particles when r > d.
 
 The two Hookean styles use this formula:
 
-.. image:: Eqs/pair_gran_hooke.jpg
-   :align: center
+.. math::
+
+   F_{hk} = (k_n \delta \mathbf{n}_{ij} -  
+   m_{eff} \gamma_n\mathbf{ v}_n) - 
+   (k_t \mathbf{ \Delta s}_t +
+   m_{eff} \gamma_t \mathbf{v}_t)
+
 
 The Hertzian style uses this formula:
 
-.. image:: Eqs/pair_gran_hertz.jpg
-   :align: center
+.. math::
+
+   F_{hz} = \sqrt{\delta} \sqrt{\frac{R_i R_j}{R_i + R_j}} F_{hk} = 
+     \sqrt{\delta} \sqrt{\frac{R_i R_j}{R_i + R_j}} 
+     \Big[ (k_n \delta \mathbf{n}_{ij} -  
+       m_{eff} \: \gamma_n \mathbf{ v}_n) -
+       (k_t \mathbf{ \Delta s}_t +
+       m_{eff} \: \gamma_t \mathbf{v}_t) \Big]
+
 
 In both equations the first parenthesized term is the normal force
 between the two particles and the second parenthesized term is the
@@ -92,34 +104,35 @@ if *dampflag* is set to 0.
 
 The other quantities in the equations are as follows:
 
-* delta = d - r = overlap distance of 2 particles
-* Kn = elastic constant for normal contact
-* Kt = elastic constant for tangential contact
-* gamma\_n = viscoelastic damping constant for normal contact
-* gamma\_t = viscoelastic damping constant for tangential contact
-* m\_eff = Mi Mj / (Mi + Mj) = effective mass of 2 particles of mass Mi and Mj
-* Delta St = tangential displacement vector between 2 particles       which is truncated to satisfy a frictional yield criterion
-* n\_ij = unit vector along the line connecting the centers of the 2 particles
-* Vn = normal component of the relative velocity of the 2 particles
-* Vt = tangential component of the relative velocity of the 2 particles
+* :math:`\delta` = d - r = overlap distance of 2 particles
+* :math:`K_n` = elastic constant for normal contact
+* :math:`K_t` = elastic constant for tangential contact
+* :math:`\gamma_n` = viscoelastic damping constant for normal contact
+* :math:`\gamma_t` = viscoelastic damping constant for tangential contact
+* :math:`m_{eff} = M_i M_j / (M_i + M_j) =` effective mass of 2 particles of mass M\_i and M\_j
+* :math:`\mathbf{\Delta s}_t =` tangential displacement vector between 2 particles       which is truncated to satisfy a frictional yield criterion
+* :math:`n_{ij} =` unit vector along the line connecting the centers of the 2 particles
+* :math:`V_n =` normal component of the relative velocity of the 2 particles
+* :math:`V_t =` tangential component of the relative velocity of the 2 particles
 
-The Kn, Kt, gamma\_n, and gamma\_t coefficients are specified as
-parameters to the pair\_style command.  If a NULL is used for Kt, then
-a default value is used where Kt = 2/7 Kn.  If a NULL is used for
-gamma\_t, then a default value is used where gamma\_t = 1/2 gamma\_n.
+The :math:`K_n`, :math:`K_t`, :math:`\gamma_n`, and :math:`\gamma_t`
+coefficients are specified as parameters to the pair\_style command.  If
+a NULL is used for :math:`K_t`, then a default value is used where
+:math:`K_t = 2/7 K_n`.  If a NULL is used for :math:`\gamma_t`, then a
+default value is used where :math:`\gamma_t = 1/2 \gamma_n`.
 
 The interpretation and units for these 4 coefficients are different in
 the Hookean versus Hertzian equations.
 
 The Hookean model is one where the normal push-back force for two
 overlapping particles is a linear function of the overlap distance.
-Thus the specified Kn is in units of (force/distance).  Note that this
-push-back force is independent of absolute particle size (in the
-monodisperse case) and of the relative sizes of the two particles (in
-the polydisperse case).  This model also applies to the other terms in
-the force equation so that the specified gamma\_n is in units of
-(1/time), Kt is in units of (force/distance), and gamma\_t is in units
-of (1/time).
+Thus the specified :math:`K_n` is in units of (force/distance).  Note
+that this push-back force is independent of absolute particle size (in
+the monodisperse case) and of the relative sizes of the two particles
+(in the polydisperse case).  This model also applies to the other terms
+in the force equation so that the specified :math:`\gamma_n` is in units
+of (1/time), :math:`K_t` is in units of (force/distance), and
+:math:`\gamma_t` is in units of (1/time).
 
 The Hertzian model is one where the normal push-back force for two
 overlapping particles is proportional to the area of overlap of the
@@ -128,38 +141,39 @@ Thus Kn has units of force per area and is thus specified in units of
 (pressure).  The effects of absolute particle size (monodispersity)
 and relative size (polydispersity) are captured in the radii-dependent
 pre-factors.  When these pre-factors are carried through to the other
-terms in the force equation it means that the specified gamma\_n is in
-units of (1/(time\*distance)), Kt is in units of (pressure), and
-gamma\_t is in units of (1/(time\*distance)).
+terms in the force equation it means that the specified :math:`\gamma_n` is in
+units of (1/(time\*distance)), :math:`K_t` is in units of (pressure), and
+:math:`\gamma_t` is in units of (1/(time\*distance)).
 
-Note that in the Hookean case, Kn can be thought of as a linear spring
-constant with units of force/distance.  In the Hertzian case, Kn is
-like a non-linear spring constant with units of force/area or
-pressure, and as shown in the :ref:`(Zhang) <Zhang3>` paper, Kn = 4G /
-(3(1-nu)) where nu = the Poisson ratio, G = shear modulus = E /
-(2(1+nu)), and E = Young's modulus.  Similarly, Kt = 4G / (2-nu).
-(NOTE: in an earlier version of the manual, we incorrectly stated that
-Kt = 8G / (2-nu).)
+Note that in the Hookean case, :math:`K_n` can be thought of as a linear
+spring constant with units of force/distance.  In the Hertzian case,
+:math:`K_n` is like a non-linear spring constant with units of
+force/area or pressure, and as shown in the :ref:`(Zhang) <Zhang3>`
+paper, :math:`K_n = 4G / (3(1-\nu))` where :math:`\nu =` the Poisson ratio,
+G = shear modulus = :math:`E / (2(1+\nu))`, and E = Young's modulus.  Similarly,
+:math:`K_t = 4G / (2-\nu)`.  (NOTE: in an earlier version of the manual, we incorrectly
+stated that :math:`K_t = 8G / (2-\nu)`.)
 
-Thus in the Hertzian case Kn and Kt can be set to values that
-corresponds to properties of the material being modeled.  This is also
-true in the Hookean case, except that a spring constant must be chosen
-that is appropriate for the absolute size of particles in the model.
-Since relative particle sizes are not accounted for, the Hookean
-styles may not be a suitable model for polydisperse systems.
+Thus in the Hertzian case :math:`K_n` and :math:`K_t` can be set to
+values that corresponds to properties of the material being modeled.
+This is also true in the Hookean case, except that a spring constant
+must be chosen that is appropriate for the absolute size of particles in
+the model.  Since relative particle sizes are not accounted for, the
+Hookean styles may not be a suitable model for polydisperse systems.
 
 .. note::
 
    In versions of LAMMPS before 9Jan09, the equation for Hertzian
-   interactions did not include the sqrt(RiRj/Ri+Rj) term and thus was
-   not as accurate for polydisperse systems.  For monodisperse systems,
-   sqrt(RiRj/Ri+Rj) is a constant factor that effectively scales all 4
-   coefficients: Kn, Kt, gamma\_n, gamma\_t.  Thus you can set the values
-   of these 4 coefficients appropriately in the current code to reproduce
+   interactions did not include the :math:`\sqrt{r_i r_j / (r_i + r_j)}`
+   term and thus was not as accurate for polydisperse systems.  For
+   monodisperse systems, :math:`\sqrt{ r_i r_j /(r_i+r_j)}` is a
+   constant factor that effectively scales all 4 coefficients:
+   :math:`K_n, K_t, \gamma_n, \gamma_t`.  Thus you can set the values of
+   these 4 coefficients appropriately in the current code to reproduce
    the results of a previous Hertzian monodisperse calculation.  For
    example, for the common case of a monodisperse system with particles
-   of diameter 1, all 4 of these coefficients should now be set 2x larger
-   than they were previously.
+   of diameter 1, all 4 of these coefficients should now be set 2x
+   larger than they were previously.
 
 Xmu is also specified in the pair\_style command and is the upper limit
 of the tangential force through the Coulomb criterion Ft = xmu\*Fn,

doc/src/pair_gromacs.rst

@@ -64,21 +64,30 @@ smoothly to zero between an inner and outer cutoff.  It is a commonly
 used potential in the `GROMACS <http://www.gromacs.org>`_ MD code and for
 the coarse-grained models of :ref:`(Marrink) <Marrink>`.
 
-.. image:: Eqs/pair_gromacs.jpg
-   :align: center
+.. math::
 
-r1 is the inner cutoff; rc is the outer cutoff.  The coefficients A, B,
-and C are computed by LAMMPS to perform the shifting and smoothing.
-The function
-S(r) is actually applied once to each term of the LJ formula and once
-to the Coulombic formula, so there are 2 or 3 sets of A,B,C coefficients
-depending on which pair\_style is used.  The boundary conditions
-applied to the smoothing function are as follows: S'(r1) = S''(r1) = 0,
-S(rc) = -E(rc), S'(rc) = -E'(rc), and S''(rc) = -E''(rc),
-where E(r) is the corresponding term
-in the LJ or Coulombic potential energy function.
-Single and double primes denote first and second
-derivatives with respect to r, respectively.
+   E_{LJ} = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - 
+            \left(\frac{\sigma}{r}\right)^6 \right] + S_{LJ}(r)
+                       \qquad r < r_c \\
+   E_C  = & \frac{C q_i q_j}{\epsilon  r} + S_C(r) \qquad r < r_c \\
+   S(r) = & C \qquad r < r_1 \\
+   S(r) = & \frac{A}{3} (r - r_1)^3 + \frac{B}{4} (r - r_1)^4 + C \qquad  r_1 < r < r_c \\
+   A = & (-3 E'(r_c) + (r_c - r_1) E''(r_c))/(r_c - r_1)^2 \\
+   B = & (2 E'(r_c) - (r_c - r_1) E''(r_c))/(r_c - r_1)^3 \\
+   C = & -E(r_c) + \frac{1}{2} (r_c - r_1) E'(r_c) - \frac{1}{12} (r_c - r_1)^2 E''(r_c)
+
+
+:math:`r_1` is the inner cutoff; :math:`r_c` is the outer cutoff.  The
+coefficients A, B, and C are computed by LAMMPS to perform the shifting
+and smoothing.  The function S(r) is actually applied once to each term
+of the LJ formula and once to the Coulombic formula, so there are 2 or 3
+sets of A,B,C coefficients depending on which pair\_style is used.  The
+boundary conditions applied to the smoothing function are as follows:
+:math:`S'(r_1) = S''(r_1) = 0, S(r_c) = -E(r_c), S'(r_c) = -E'(r_c)`,
+and :math:`S''(r_c) = -E''(r_c)`, where E(r) is the corresponding term
+in the LJ or Coulombic potential energy function.  Single and double
+primes denote first and second derivatives with respect to r,
+respectively.
 
 The inner and outer cutoff for the LJ and Coulombic terms can be the
 same or different depending on whether 2 or 4 arguments are used in
@@ -91,14 +100,13 @@ above, or in the data file or restart files read by the
 :doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
 commands, or by mixing as described below:
 
-* epsilon (energy units)
-* sigma (distance units)
+* :math:`\epsilon` (energy units)
+* :math:`\sigma` (distance units)
 * inner (distance units)
 * outer (distance units)
 
 Note that sigma is defined in the LJ formula as the zero-crossing
-distance for the potential, not as the energy minimum at 2\^(1/6)
-sigma.
+distance for the potential, not as the energy minimum at :math:`2^{1/6} \sigma`.
 
 The last 2 coefficients are optional inner and outer cutoffs for style
 *lj/gromacs*\ .  If not specified, the global *inner* and *outer* values