Shell molecular dynamics
------------------------
|Gromacs| can simulate polarizability using the shell model of Dick and
OverhauserĀ \ :ref:`43 `. In such models a shell particle
representing the electronic degrees of freedom is attached to a nucleus
by a spring. The potential energy is minimized with respect to the shell
position at every step of the simulation (see below). Successful
applications of shell models in |Gromacs| have been published for
:math:`N_2` :ref:`44 ` and water\ :ref:`45 `.
Optimization of the shell positions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The force :math:`\mathbf{F}`\ :math:`_S` on a shell
particle :math:`S` can be decomposed into two components
.. math:: \mathbf{F}_S ~=~ \mathbf{F}_{bond} + \mathbf{F}_{nb}
:label: eqnshellforcedecomp
where :math:`\mathbf{F}_{bond}` denotes the
component representing the polarization energy, usually represented by a
harmonic potential and :math:`\mathbf{F}_{nb}` is the sum of Coulomb
and van der Waals interactions. If we assume that
:math:`\mathbf{F}_{nb}` is almost constant we
can analytically derive the optimal position of the shell, i.e. where
:math:`\mathbf{F}_S` = 0. If we have the shell S connected to atom A we have
.. math:: \mathbf{F}_{bond} ~=~ k_b \left( \mathbf{x}_S - \mathbf{x}_A\right).
:label: eqnshell
In an iterative solver, we have positions :math:`\mathbf{x}_S(n)` where :math:`n` is
the iteration count. We now have at iteration :math:`n`
.. math:: \mathbf{F}_{nb} ~=~ \mathbf{F}_S - k_b \left( \mathbf{x}_S(n) - \mathbf{x}_A\right)
:label: eqnshellsolv
and the optimal position for the shells :math:`x_S(n+1)` thus follows from
.. math:: \mathbf{F}_S - k_b \left( \mathbf{x}_S(n) - \mathbf{x}_A\right) + k_b \left( \mathbf{x}_S(n+1) - \mathbf{x}_A\right) = 0
:label: eqnshelloptpos
if we write
.. math:: \Delta \mathbf{x}_S = \mathbf{x}_S(n+1) - \mathbf{x}_S(n)
:label: eqnshelloptpos2
we finally obtain
.. math:: \Delta \mathbf{x}_S = \mathbf{F}_S/k_b
:label: eqnshelloptpos3
which then yields the algorithm to compute the next trial in the
optimization of shell positions
.. math:: \mathbf{x}_S(n+1) ~=~ \mathbf{x}_S(n) + \mathbf{F}_S/k_b.
:label: eqnshelloptpos4