remove \begin{equation} \end{equation} which are not needed and break epub

have to be modified to include the Drude dipoles and how to handle

them.

----------

**Overview of Drude induced dipoles**

Polarizable atoms acquire an induced electric dipole moment under the

.. math::

   \begin{equation} K_D = \frac 1 2\, \frac {q_D^2} \alpha\end{equation}

   K_D = \frac 1 2\, \frac {q_D^2} \alpha

Ideally, the mass of the Drude particle should be small, and the

stiffness of the harmonic bond should be large, so that the Drude

#. The possibility to thermostat the additional degrees of freedom associated with the induced dipoles at very low temperature, in terms of the reduced coordinates of the Drude particles with respect to their cores. This makes the trajectory close to that of relaxed induced dipoles.

#. The Drude dipoles on covalently bonded atoms interact too strongly due to the short distances, so an atom may capture the Drude particle (shell) of a neighbor, or the induced dipoles within the same molecule may align too much.  To avoid this, damping at short of the interactions between the point charges composing the induced dipole can be done by :ref:`Thole <Thole2>` functions.

----------

**Preparation of the data file**

The data file is similar to a standard LAMMPS data file for

.. math::

   \begin{equation}\left< U\right>_{HMA} = \frac{d}{2} (N-1) k_B T  + \left< U + \frac{1}{2} F\bullet\Delta r \right>\end{equation}

   \left< U\right>_{HMA} = \frac{d}{2} (N-1) k_B T  + \left< U + \frac{1}{2} F\bullet\Delta r \right>

where :math:`N` is the number of atoms in the system, :math:`k_B` is Boltzmann's

constant, :math:`T` is the temperature, :math:`d` is the

.. math::

   \begin{equation}\left< P\right>_{HMA} = \Delta \hat P + \left< P_{vir} + \frac{\beta \Delta \hat P - \rho}{d(N-1)} F\bullet\Delta r \right>\end{equation}

   \left< P\right>_{HMA} = \Delta \hat P + \left< P_{vir} + \frac{\beta \Delta \hat P - \rho}{d(N-1)} F\bullet\Delta r \right>

where :math:`\rho` is the number density of the system, :math:`\Delta \hat P` is the

difference between the harmonic and lattice pressure, :math:`P_{vir}` is

.. math::

   \begin{equation}\left<C_V \right>_{HMA} = \frac{d}{2} (N-1) k_B + \frac{1}{k_B T^2} \left( \left<

   \left<C_V \right>_{HMA} = \frac{d}{2} (N-1) k_B + \frac{1}{k_B T^2} \left( \left<

   U_{HMA}^2 \right> - \left<U_{HMA}\right>^2 \right) + \frac{1}{4 T}

   \left< F\bullet\Delta r + \Delta r \bullet \Phi \bullet \Delta r \right>\end{equation}

   \left< F\bullet\Delta r + \Delta r \bullet \Phi \bullet \Delta r \right>

where :math:`\Phi` is the Hessian matrix. The compute hma command

computes the full expression for :math:`C_V` except for the

.. math::

   \begin{equation} M' = M + m \end{equation}

    M' = M + m 

.. math::

   \begin{equation} m' = \frac {M\, m } {M'} \end{equation}

    m' = \frac {M\, m } {M'} 

Positions:

.. math::

   \begin{equation} X' = \frac {M\, X + m\, x} {M'}\end{equation}

    X' = \frac {M\, X + m\, x} {M'}

.. math::

   \begin{equation} x' = x - X \end{equation}

    x' = x - X 

Velocities:

.. math::

   \begin{equation} V' = \frac {M\, V + m\, v} {M'}\end{equation}

    V' = \frac {M\, V + m\, v} {M'}

.. math::

   \begin{equation} v' = v - V \end{equation}

    v' = v - V 

Forces:

.. math::

   \begin{equation} F' = F + f \end{equation}

    F' = F + f 

.. math::

   \begin{equation} f' = \frac { M\, f - m\, F} {M'}\end{equation}

    f' = \frac { M\, f - m\, F} {M'}

This transform conserves the total kinetic energy

.. math::

   \begin{equation} \frac 1 2 \, (M\, V^2\ + m\, v^2)

   = \frac 1 2 \, (M'\, V'^2\ + m'\, v'^2) \end{equation}

    \frac 1 2 \, (M\, V^2\ + m\, v^2)

   = \frac 1 2 \, (M'\, V'^2\ + m'\, v'^2) 

and the virial defined with absolute positions

.. math::

   \begin{equation} X\, F + x\, f = X'\, F' + x'\, f' \end{equation}

    X\, F + x\, f = X'\, F' + x'\, f' 

----------

.. math::

   \begin{equation}\vec{F}_i = \vec{F}^0_i - \frac{\vec{v}_i}{\|\vec{v}_i\|} \cdot S_e\end{equation}

   \vec{F}_i = \vec{F}^0_i - \frac{\vec{v}_i}{\|\vec{v}_i\|} \cdot S_e

where :math:`\vec{F}_i` is the resulting total force on the atom.

:math:`\vec{F}^0_i` is the original force applied to the atom, :math:`\vec{v}_i` is

.. math::

   \begin{equation} \frac {dq}{dt} = \frac{p}{m}, \end{equation}

    \frac {dq}{dt} = \frac{p}{m},

.. math::

   \begin{equation} \frac {dp}{dt} = -\gamma p + W + F, \end{equation}

   \frac {dp}{dt} = -\gamma p + W + F,

where :math:`F` is the physical force, :math:`\gamma` is the friction coefficient, and :math:`W` is a

Gaussian random force.