first chunk of compute commands to be converted to use embedded math

doc/src/Eqs/centro_symmetry.tex

@@ -1,9 +0,0 @@
-\documentclass[12pt]{article}
-
-\begin{document}
-
-$$
-   CS = \sum_{i = 1}^{N/2} | \vec{R}_i + \vec{R}_{i+N/2} |^2
-$$
-
-\end{document}

doc/src/Eqs/cna_cutoff1.tex

@@ -1,14 +0,0 @@
-\documentclass[12pt,article]{article}
-
-\usepackage{indentfirst}
-\usepackage{amsmath}
-
-\begin{document}
-
-\begin{eqnarray*}
-  r_{c}^{fcc} & = & \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) \mathrm{a} \simeq 0.8536 \:\mathrm{a} \\
-  r_{c}^{bcc} & = & \frac{1}{2}(\sqrt{2} + 1) \mathrm{a} \simeq 1.207 \:\mathrm{a} \\
-  r_{c}^{hcp} & = & \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) \mathrm{a}
-\end{eqnarray*}
-
-\end{document}

doc/src/Eqs/cna_cutoff2.tex

@@ -1,12 +0,0 @@
-\documentclass[12pt,article]{article}
-
-\usepackage{indentfirst}
-\usepackage{amsmath}
-
-\begin{document}
-
-$$
-  Rc + Rs > 2*{\rm cutoff}
-$$
-
-\end{document}

doc/src/Eqs/cnp_cutoff2.tex

@@ -1,12 +0,0 @@
-\documentclass[12pt,article]{article}
-
-\usepackage{indentfirst}
-\usepackage{amsmath}
-
-\begin{document}
-
-$$
-  Rc + Rs > 2*{\rm cutoff}
-$$
-
-\end{document}

doc/src/Eqs/cnp_eq.tex

@@ -1,9 +0,0 @@
-\documentclass[12pt]{article}
-
-\begin{document}
-
-$$
-   Q_{i} = \frac{1}{n_i}\sum_{j = 1}^{n_i} | \sum_{k = 1}^{n_{ij}}  \vec{R}_{ik} + \vec{R}_{jk} |^2
-$$
-
-\end{document}

doc/src/Eqs/compute_gyration.tex

@@ -1,9 +0,0 @@
-\documentstyle[12pt]{article}
-
-\begin{document}
-
-$$
- {R_g}^2 = \frac{1}{M} \sum_i m_i (r_i - r_{cm})^2
-$$
-
-\end{document}

doc/src/Eqs/compute_shape_parameters.tex

@@ -1,13 +0,0 @@
-\documentclass[12pt]{article}
-
-\pagestyle{empty}
-\begin{document}
-
-\begin{eqnarray*}
- c = l_z - 0.5(l_y+l_x) \\
- b = l_y - l_x \\
- k = \frac{3}{2} \frac{l_x^2+l_y^2+l_z^2}{(l_x+l_y+l_z)^2} - \frac{1}{2} 
-\end{eqnarray*}
-
-\end{document}
-

doc/src/compute_centro_atom.rst

@@ -52,20 +52,23 @@ in the specified compute group.
 This parameter is computed using the following formula from
 :ref:`(Kelchner) <Kelchner>`
 
-.. image:: Eqs/centro_symmetry.jpg
-   :align: center
+.. math::
 
-where the *N* nearest neighbors of each atom are identified and Ri and
-Ri+N/2 are vectors from the central atom to a particular pair of
-nearest neighbors.  There are N\*(N-1)/2 possible neighbor pairs that
-can contribute to this formula.  The quantity in the sum is computed
-for each, and the N/2 smallest are used.  This will typically be for
-pairs of atoms in symmetrically opposite positions with respect to the
-central atom; hence the i+N/2 notation.
+   CS = \sum_{i = 1}^{N/2} | \vec{R}_i + \vec{R}_{i+N/2} |^2
 
-*N* is an input parameter, which should be set to correspond to the
-number of nearest neighbors in the underlying lattice of atoms.  If
-the keyword *fcc* or *bcc* is used, *N* is set to 12 and 8
+
+where the :math:`N` nearest neighbors of each atom are identified and
+:math:`\vec{R}_i` and :math:`\vec{R}_{i+N/2}` are vectors from the
+central atom to a particular pair of nearest neighbors.  There are
+:math:`N (N-1)/2` possible neighbor pairs that can contribute to this
+formula.  The quantity in the sum is computed for each, and the
+:math:`N/2` smallest are used.  This will typically be for pairs of
+atoms in symmetrically opposite positions with respect to the central
+atom; hence the :math:`i+N/2` notation.
+
+:math:`N` is an input parameter, which should be set to correspond to
+the number of nearest neighbors in the underlying lattice of atoms.
+If the keyword *fcc* or *bcc* is used, *N* is set to 12 and 8
 respectively.  More generally, *N* can be set to a positive, even
 integer.
 
@@ -74,9 +77,9 @@ lattice, the centro-symmetry parameter will be 0.  It will be near 0
 for small thermal perturbations of a perfect lattice.  If a point
 defect exists, the symmetry is broken, and the parameter will be a
 larger positive value.  An atom at a surface will have a large
-positive parameter.  If the atom does not have *N* neighbors (within
-the potential cutoff), then its centro-symmetry parameter is set to
-0.0.
+positive parameter.  If the atom does not have :math:`N` neighbors
+(within the potential cutoff), then its centro-symmetry parameter is
+set to 0.0.
 
 If the keyword *axes* has the setting *yes*\ , then this compute also
 estimates three symmetry axes for each atom's local neighborhood.  The
@@ -95,7 +98,7 @@ of any atom.
 
 Only atoms within the cutoff of the pairwise neighbor list are
 considered as possible neighbors.  Atoms not in the compute group are
-included in the *N* neighbors used in this calculation.
+included in the :math:`N` neighbors used in this calculation.
 
 The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (e.g. each time a snapshot of atoms

doc/src/compute_cna_atom.rst

@@ -51,8 +51,12 @@ E.g. 12 nearest neighbor for perfect FCC and HCP crystals, 14 nearest
 neighbors for perfect BCC crystals.  These formulas can be used to
 obtain a good cutoff distance:
 
-.. image:: Eqs/cna_cutoff1.jpg
-   :align: center
+.. math::
+
+  r_{c}^{fcc} = & \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) \mathrm{a} \simeq 0.8536 \:\mathrm{a} \\
+  r_{c}^{bcc} = & \frac{1}{2}(\sqrt{2} + 1) \mathrm{a} \simeq 1.207 \:\mathrm{a} \\
+  r_{c}^{hcp} = & \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) \mathrm{a}
+
 
 where a is the lattice constant for the crystal structure concerned
 and in the HCP case, x = (c/a) / 1.633, where 1.633 is the ideal c/a
@@ -62,10 +66,13 @@ Also note that since the CNA calculation in LAMMPS uses the neighbors
 of an owned atom to find the nearest neighbors of a ghost atom, the
 following relation should also be satisfied:
 
-.. image:: Eqs/cna_cutoff2.jpg
-   :align: center
+.. math::
 
-where Rc is the cutoff distance of the potential, Rs is the skin
+  r_c + r_s > 2*{\rm cutoff}
+
+
+where :math:`r_c` is the cutoff distance of the potential, :math:`r_s`
+is the skin
 distance as specified by the :doc:`neighbor <neighbor>` command, and
 cutoff is the argument used with the compute cna/atom command.  LAMMPS
 will issue a warning if this is not the case.

doc/src/compute_cnp_atom.rst

@@ -40,13 +40,16 @@ only be performed on single component systems.
 This parameter is computed using the following formula from
 :ref:`(Tsuzuki) <Tsuzuki2>`
 
-.. image:: Eqs/cnp_eq.jpg
-   :align: center
+.. math::
 
-where the index *j* goes over the *n*\ i nearest neighbors of atom
-*i*\ , and the index *k* goes over the *n*\ ij common nearest neighbors
-between atom *i* and atom *j*\ . Rik and Rjk are the vectors connecting atom
-*k* to atoms *i* and *j*\ .  The quantity in the double sum is computed
+   Q_{i} = \frac{1}{n_i}\sum_{j = 1}^{n_i} | \sum_{k = 1}^{n_{ij}}  \vec{R}_{ik} + \vec{R}_{jk} |^2
+
+
+where the index *j* goes over the :math:`n_i` nearest neighbors of atom
+*i*\ , and the index *k* goes over the :math:`n_{ij}` common nearest neighbors
+between atom *i* and atom *j*\ . :math:`\vec{R}_{ik}` and
+:math:`\vec{R}_{jk}` are the vectors connecting atom *k* to atoms *i*
+and *j*\ .  The quantity in the double sum is computed
 for each atom.
 
 The CNP calculation is sensitive to the specified cutoff value.
@@ -56,8 +59,12 @@ E.g. 12 nearest neighbor for perfect FCC and HCP crystals, 14 nearest
 neighbors for perfect BCC crystals.  These formulas can be used to
 obtain a good cutoff distance:
 
-.. image:: Eqs/cnp_cutoff.jpg
-   :align: center
+.. math::
+
+  r_{c}^{fcc} = & \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) \mathrm{a} \simeq 0.8536 \:\mathrm{a} \\
+  r_{c}^{bcc} = & \frac{1}{2}(\sqrt{2} + 1) \mathrm{a} \simeq 1.207 \:\mathrm{a} \\
+  r_{c}^{hcp} = & \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) \mathrm{a}
+
 
 where a is the lattice constant for the crystal structure concerned
 and in the HCP case, x = (c/a) / 1.633, where 1.633 is the ideal c/a
@@ -67,10 +74,13 @@ Also note that since the CNP calculation in LAMMPS uses the neighbors
 of an owned atom to find the nearest neighbors of a ghost atom, the
 following relation should also be satisfied:
 
-.. image:: Eqs/cnp_cutoff2.jpg
-   :align: center
+.. math::
 
-where Rc is the cutoff distance of the potential, Rs is the skin
+  r_c + r_s > 2*{\rm cutoff}
+
+
+where :math:`r_c` is the cutoff distance of the potential, :math:`r_s` is
+the skin
 distance as specified by the :doc:`neighbor <neighbor>` command, and
 cutoff is the argument used with the compute cnp/atom command.  LAMMPS
 will issue a warning if this is not the case.

doc/src/compute_dpd.rst

@@ -26,19 +26,26 @@ Description
 """""""""""
 
 Define a computation that accumulates the total internal conductive
-energy (U\_cond), the total internal mechanical energy (U\_mech), the
-total chemical energy (U\_chem) and the *harmonic* average of the internal
-temperature (dpdTheta) for the entire system of particles.  See the
+energy (:math:`U^{cond}`), the total internal mechanical energy
+(:math:`U^{mech}`), the total chemical energy (:math:`U^{chem}`)
+and the *harmonic* average of the internal temperature (:math:`\theta_{avg}`)
+for the entire system of particles.  See the
 :doc:`compute dpd/atom <compute_dpd_atom>` command if you want
 per-particle internal energies and internal temperatures.
 
 The system internal properties are computed according to the following
 relations:
 
-.. image:: Eqs/compute_dpd.jpg
-   :align: center
+.. math::
 
-where N is the number of particles in the system
+   U^{cond} = & \displaystyle\sum_{i=1}^{N} u_{i}^{cond} \\ 
+   U^{mech} = & \displaystyle\sum_{i=1}^{N} u_{i}^{mech} \\
+   U^{chem} = & \displaystyle\sum_{i=1}^{N} u_{i}^{chem} \\
+          U = & \displaystyle\sum_{i=1}^{N} (u_{i}^{cond} + u_{i}^{mech} + u_{i}^{chem}) \\
+   \theta_{avg} = & (\frac{1}{N}\displaystyle\sum_{i=1}^{N} \frac{1}{\theta_{i}})^{-1} \\
+
+
+where :math:`N` is the number of particles in the system
 
 
 ----------
@@ -46,8 +53,9 @@ where N is the number of particles in the system
 
 **Output info:**
 
-This compute calculates a global vector of length 5 (U\_cond, U\_mech,
-U\_chem, dpdTheta, N\_particles), which can be accessed by indices 1-5.
+This compute calculates a global vector of length 5 (:math:`U^{cond}`,
+:math:`U^{mech}`, :math:`U^{chem}`, :math:`\theta_{avg}`, :math:`N`),
+which can be accessed by indices 1-5.
 See the :doc:`Howto output <Howto_output>` doc page for an overview of
 LAMMPS output options.
 

doc/src/compute_dpd_atom.rst

@@ -23,10 +23,10 @@ Description
 """""""""""
 
 Define a computation that accesses the per-particle internal
-conductive energy (u\_cond), internal mechanical energy (u\_mech),
-internal chemical energy (u\_chem) and
-internal temperatures (dpdTheta) for each particle in a group.  See
-the :doc:`compute dpd <compute_dpd>` command if you want the total
+conductive energy (:math:`u^{cond}`), internal mechanical
+energy (:math:`u^{mech}`), internal chemical energy (:math:`u^{chem}`)
+and internal temperatures (:math:`\theta`) for each particle in a group.
+See the :doc:`compute dpd <compute_dpd>` command if you want the total
 internal conductive energy, the total internal mechanical energy, the
 total chemical energy and
 average internal temperature of the entire system or group of dpd
@@ -34,14 +34,16 @@ particles.
 
 **Output info:**
 
-This compute calculates a per-particle array with 4 columns (u\_cond,
-u\_mech, u\_chem, dpdTheta), which can be accessed by indices 1-4 by any
+This compute calculates a per-particle array with 4 columns (:math:`u^{cond}`,
+:math:`u^{mech}`, :math:`u^{chem}`, :math:`\theta`), which can be accessed
+by indices 1-4 by any
 command that uses per-particle values from a compute as input.  See
 the :doc:`Howto output <Howto_output>` doc page for an overview of
 LAMMPS output options.
 
-The per-particle array values will be in energy (u\_cond, u\_mech, u\_chem)
-and temperature (dpdTheta) :doc:`units <units>`.
+The per-particle array values will be in energy (:math:`u^{cond}`,
+:math:`u^{mech}`, :math:`u^{chem}`)
+and temperature (:math:`theta`) :doc:`units <units>`.
 
 Restrictions
 """"""""""""

doc/src/compute_entropy_atom.rst

@@ -53,27 +53,33 @@ information about the solid structure is required.
 This parameter for atom i is computed using the following formula from
 :ref:`(Piaggi) <Piaggi>` and :ref:`(Nettleton) <Nettleton>` ,
 
-.. image:: Eqs/pair_entropy.jpg
-   :align: center
+.. math::
+
+   s_S^i=-2\pi\rho k_B \int\limits_0^{r_m} \left [ g(r) \ln g(r) - g(r) + 1 \right ] r^2 dr
+
 
 where r is a distance, g(r) is the radial distribution function of atom
 i and rho is the density of the system. The g(r) computed for each
 atom i can be noisy and therefore it is smoothed using:
 
-.. image:: Eqs/pair_entropy2.jpg
-   :align: center
+.. math::
 
-where the sum in j goes through the neighbors of atom i, and sigma is a
-parameter to control the smoothing.
+   g_m^i(r) = \frac{1}{4 \pi \rho r^2} \sum\limits_{j} \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(r-r_{ij})^2/(2\sigma^2)}
 
-The input parameters are *sigma* the smoothing parameter, and the
-*cutoff* for the calculation of g(r).
+
+where the sum in j goes through the neighbors of atom i, and :math:`\sigma`
+is a parameter to control the smoothing.
+
+The input parameters are *sigma* the smoothing parameter :math:`\sigma`,
+and the *cutoff* for the calculation of g(r).
 
 If the keyword *avg* has the setting *yes*\ , then this compute also
 averages the parameter over the neighbors  of atom i according to:
 
-.. image:: Eqs/pair_entropy3.jpg
-   :align: center
+.. math::
+
+  \left< s_S^i \right>  = \frac{\sum_j s_S^j + s_S^i}{N + 1}
+
 
 where the sum j goes over the neighbors of atom i and N is the number
 of neighbors. This procedure provides a sharper distinction between

doc/src/compute_fep.rst

@@ -70,14 +70,18 @@ initial interactions of the atoms that will undergo perturbation, and
 a term :math:`U_1` corresponding to the final interactions of
 these atoms:
 
-.. image:: Eqs/compute_fep_u.jpg
-   :align: center
+.. math::
+
+   U(\lambda) = U_{\mathrm{bg}} + U_1(\lambda) + U_0(\lambda)
 
 A coupling parameter :math:`\lambda` varying from 0 to 1 connects the
 reference and perturbed systems:
 
-.. image:: Eqs/compute_fep_lambda.jpg
-   :align: center
+.. math::
+
+  \lambda = 0 \quad\Rightarrow\quad U = U_{\mathrm{bg}} + U_0 \\
+  \lambda = 1 \quad\Rightarrow\quad U = U_{\mathrm{bg}} + U_1 
+
 
 It is possible but not necessary that the coupling parameter (or a
 function thereof) appears as a multiplication factor of the potential
@@ -89,16 +93,23 @@ This command can be combined with :doc:`fix adapt <fix_adapt>` to
 perform multistage free-energy perturbation calculations along
 stepwise alchemical transformations during a simulation run:
 
-.. image:: Eqs/compute_fep_fep.jpg
-   :align: center
+.. math::
+
+   \Delta_0^1 A = \sum_{i=0}^{n-1} \Delta_{\lambda_i}^{\lambda_{i+1}} A =
+   - kT \sum_{i=0}^{n-1} \ln \left< \exp \left( - \frac{U(\lambda_{i+1}) -
+      U(\lambda_i)}{kT} \right) \right>_{\lambda_i}
 
 This compute is suitable for the finite-difference thermodynamic
 integration (FDTI) method :ref:`(Mezei) <Mezei>`, which is based on an
 evaluation of the numerical derivative of the free energy by a
 perturbation method using a very small :math:`\delta`:
 
-.. image:: Eqs/compute_fep_fdti.jpg
-   :align: center
+.. math::
+
+   \Delta_0^1 A = \int_{\lambda=0}^{\lambda=1} \left( \frac{\partial
+   A(\lambda)}{\partial\lambda} \right)_\lambda \mathrm{d}\lambda
+   \approx \sum_{i=0}^{n-1} w_i \frac{A(\lambda_{i} + \delta) -
+   A(\lambda_i)}{\delta}
 
 where :math:`w_i` are weights of a numerical quadrature. The :doc:`fix adapt <fix_adapt>` command can be used to define the stages of
 :math:`\lambda` at which the derivative is calculated and averaged.
@@ -109,16 +120,23 @@ choosing a very small perturbation :math:`\delta` the thermodynamic
 integration method can be implemented using a numerical evaluation of
 the derivative of the potential energy with respect to :math:`\lambda`:
 
-.. image:: Eqs/compute_fep_ti.jpg
-   :align: center
+.. math::
+
+   \Delta_0^1 A = \int_{\lambda=0}^{\lambda=1} \left< \frac{\partial
+   U(\lambda)}{\partial\lambda} \right>_\lambda \mathrm{d}\lambda
+   \approx \sum_{i=0}^{n-1} w_i \left< \frac{U(\lambda_{i} + \delta) -
+   U(\lambda_i)}{\delta} \right>_{\lambda_i}
+
+
 
 Another technique to calculate free energy differences is the
 acceptance ratio method :ref:`(Bennet) <Bennet>`, which can be implemented
 by calculating the potential energy differences with :math:`\delta` = 1.0 on
 both the forward and reverse routes:
 
-.. image:: Eqs/compute_fep_bar.jpg
-   :align: center
+.. math::
+
+   \left< \frac{1}{1 + \exp\left[\left(U_1 - U_0 - \Delta_0^1A \right) /kT \right]} \right>_0 = \left< \frac{1}{1 + \exp\left[\left(U_0 - U_1 + \Delta_0^1A \right) /kT \right]} \right>_1
 
 The value of the free energy difference is determined by numerical
 root finding to establish the equality.
@@ -265,9 +283,11 @@ If the keyword *volume* = *yes*\ , then the Boltzmann term is multiplied
 by the volume so that correct ensemble averaging can be performed over
 trajectories during which the volume fluctuates or changes :ref:`(Allen and Tildesley) <AllenTildesley>`:
 
-.. image:: Eqs/compute_fep_vol.jpg
-   :align: center
+.. math::
 
+   \Delta_0^1 A = - kT \sum_{i=0}^{n-1} \ln \frac{\left< V \exp \left( -
+   \frac{U(\lambda_{i+1}) - U(\lambda_i)}{kT} \right)
+    \right>_{\lambda_i}}{\left< V \right>_{\lambda_i}}
 
 ----------
 

doc/src/compute_gyration.rst

@@ -32,24 +32,27 @@ periodic boundaries.
 Rg is a measure of the size of the group of atoms, and is computed as
 the square root of the Rg\^2 value in this formula
 
-.. image:: Eqs/compute_gyration.jpg
-   :align: center
+.. math::
 
-where M is the total mass of the group, Rcm is the center-of-mass
-position of the group, and the sum is over all atoms in the group.
+ {R_g}^2 = \frac{1}{M} \sum_i m_i (r_i - r_{cm})^2
 
-A Rg\^2 tensor, stored as a 6-element vector, is also calculated by
-this compute.  The formula for the components of the tensor is the
-same as the above formula, except that (Ri - Rcm)\^2 is replaced by
-(Rix - Rcmx) \* (Riy - Rcmy) for the xy component, etc.  The 6
-components of the vector are ordered xx, yy, zz, xy, xz, yz.  Note
-that unlike the scalar Rg, each of the 6 values of the tensor is
-effectively a "squared" value, since the cross-terms may be negative
+
+where :math:`M` is the total mass of the group, :math:`r_{cm}` is the
+center-of-mass position of the group, and the sum is over all atoms in
+the group.
+
+A :math:`{R_g}^2` tensor, stored as a 6-element vector, is also calculated
+by this compute.  The formula for the components of the tensor is the
+same as the above formula, except that :math:`(r_i - r_{cm})^2` is replaced
+by :math:`(r_{i,x} - r_{cm,x}) \cdot (r_{i,y} - r_{cm,y})` for the xy component,
+and so on.  The 6 components of the vector are ordered xx, yy, zz, xy, xz, yz.
+Note that unlike the scalar :math:`R_g`, each of the 6 values of the tensor
+is effectively a "squared" value, since the cross-terms may be negative
 and taking a sqrt() would be invalid.
 
 .. note::
 
-   The coordinates of an atom contribute to Rg in "unwrapped" form,
+   The coordinates of an atom contribute to :math:`R_g` in "unwrapped" form,
    by using the image flags associated with each atom.  See the :doc:`dump custom <dump>` command for a discussion of "unwrapped" coordinates.
    See the Atoms section of the :doc:`read_data <read_data>` command for a
    discussion of image flags and how they are set for each atom.  You can
@@ -58,8 +61,8 @@ and taking a sqrt() would be invalid.
 
 **Output info:**
 
-This compute calculates a global scalar (Rg) and a global vector of
-length 6 (Rg\^2 tensor), which can be accessed by indices 1-6.  These
+This compute calculates a global scalar (:math:`R_g`) and a global vector of
+length 6 (:math:`{R_g}^2` tensor), which can be accessed by indices 1-6.  These
 values can be used by any command that uses a global scalar value or
 vector values from a compute as input.  See the :doc:`Howto output <Howto_output>` doc page for an overview of LAMMPS output
 options.

doc/src/compute_gyration_chunk.rst

@@ -52,11 +52,13 @@ boundaries.
 Rg is a measure of the size of a chunk, and is computed by this
 formula
 
-.. image:: Eqs/compute_gyration.jpg
-   :align: center
+.. math::
 
-where M is the total mass of the chunk, Rcm is the center-of-mass
-position of the chunk, and the sum is over all atoms in the
+ {R_g}^2 = \frac{1}{M} \sum_i m_i (r_i - r_{cm})^2
+
+
+where :math:`M` is the total mass of the chunk, :math:`r_{cm}` is
+the center-of-mass position of the chunk, and the sum is over all atoms in the
 chunk.
 
 Note that only atoms in the specified group contribute to the
@@ -70,14 +72,16 @@ non-zero chunk IDs.
 If the *tensor* keyword is specified, then the scalar Rg value is not
 calculated, but an Rg tensor is instead calculated for each chunk.
 The formula for the components of the tensor is the same as the above
-formula, except that (Ri - Rcm)\^2 is replaced by (Rix - Rcmx) \* (Riy -
-Rcmy) for the xy component, etc.  The 6 components of the tensor are
+formula, except that :math:`(r_i - r_{cm})^2` is replaced by
+:math:`(r_{i,x} - r_{cm,x}) \cdot (r_{i,y} - r_{cm,y})` for the xy
+component, and so on.  The 6 components of the tensor are
 ordered xx, yy, zz, xy, xz, yz.
 
 .. note::
 
-   The coordinates of an atom contribute to Rg in "unwrapped" form,
-   by using the image flags associated with each atom.  See the :doc:`dump custom <dump>` command for a discussion of "unwrapped" coordinates.
+   The coordinates of an atom contribute to :math:`R_g` in "unwrapped" form,
+   by using the image flags associated with each atom.  See the :doc:`dump custom <dump>`
+   command for a discussion of "unwrapped" coordinates.
    See the Atoms section of the :doc:`read_data <read_data>` command for a
    discussion of image flags and how they are set for each atom.  You can
    reset the image flags (e.g. to 0) before invoking this compute by

doc/src/compute_gyration_shape.rst

@@ -33,10 +33,14 @@ due to atoms passing through periodic boundaries.
 The three computed shape parameters are the asphericity, b, the acylindricity, c,
 and the relative shape anisotropy, k:
 
-.. image:: Eqs/compute_shape_parameters.jpg
-   :align: center
+.. math::
 
-where lx <= ly <= lz are the three eigenvalues of the gyration tensor. A general description 
+ c = & l_z - 0.5(l_y+l_x) \\
+ b = & l_y - l_x \\
+ k = & \frac{3}{2} \frac{l_x^2+l_y^2+l_z^2}{(l_x+l_y+l_z)^2} - \frac{1}{2} 
+
+
+where :math:`l_x` <= :math:`l_y` <= :math:`l_z` are the three eigenvalues of the gyration tensor. A general description 
 of these parameters is provided in :ref:`(Mattice) <Mattice1>` while an application to polymer systems 
 can be found in :ref:`(Theodorou) <Theodorou1>`.
 The asphericity  is always non-negative and zero only when the three principal

doc/src/compute_gyration_shape_chunk.rst

@@ -33,10 +33,14 @@ all effects due to atoms passing through periodic boundaries.
 The three computed shape parameters are the asphericity, b, the acylindricity, c,
 and the relative shape anisotropy, k:
 
-.. image:: Eqs/compute_shape_parameters.jpg
-   :align: center
+.. math::
 
-where lx <= ly <= lz are the three eigenvalues of the gyration tensor. A general description 
+ c = & l_z - 0.5(l_y+l_x) \\
+ b = & l_y - l_x \\
+ k = & \frac{3}{2} \frac{l_x^2+l_y^2+l_z^2}{(l_x+l_y+l_z)^2} - \frac{1}{2} 
+
+
+where :math:`l_x` <= :math:`l_y` <= :math`l_z` are the three eigenvalues of the gyration tensor. A general description 
 of these parameters is provided in :ref:`(Mattice) <Mattice2>` while an application to polymer systems 
 can be found in :ref:`(Theodorou) <Theodorou2>`. The asphericity  is always non-negative and zero 
 only when the three principal moments are equal. This zero condition is met when the distribution