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}_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