docs/reference-manual/algorithms/shell-molecular-dynamics.rst

   1 Shell molecular dynamics
   2 ------------------------
   3
   4 |Gromacs| can simulate polarizability using the shell model of Dick and
   5 Overhauser \ :ref:`43 <refDick58>`. In such models a shell particle
   6 representing the electronic degrees of freedom is attached to a nucleus
   7 by a spring. The potential energy is minimized with respect to the shell
   8 position at every step of the simulation (see below). Successful
   9 applications of shell models in |Gromacs| have been published for
  10 :math:`N_2` :ref:`44 <refJordan95>` and water\ :ref:`45 <refMaaren2001a>`.
  11
  12 Optimization of the shell positions
  13 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  14
  15 The force
  16 :math:`\mathbf{F}`\ :math:`_S` on a shell
  17 particle :math:`S` can be decomposed into two components
  18
  19 .. math:: \mathbf{F}_S ~=~ \mathbf{F}_{bond} + \mathbf{F}_{nb}
  20
  21 where
  22 :math:`\mathbf{F}_{bond}` denotes the
  23 component representing the polarization energy, usually represented by a
  24 harmonic potential and
  25 :math:`\mathbf{F}_{nb}` is the sum of Coulomb
  26 and van der Waals interactions. If we assume that
  27 :math:`\mathbf{F}_{nb}` is almost constant we
  28 can analytically derive the optimal position of the shell, i.e. where
  29 :math:`\mathbf{F}_S` = 0. If we have the
  30 shell S connected to atom A we have
  31
  32 .. math:: \mathbf{F}_{bond} ~=~ k_b \left( \mathbf{x}_S - \mathbf{x}_A\right).
  33
  34 In an iterative solver, we have positions
  35 :math:`\mathbf{x}_S(n)` where :math:`n` is
  36 the iteration count. We now have at iteration :math:`n`
  37
  38 .. math:: \mathbf{F}_{nb} ~=~ \mathbf{F}_S - k_b \left( \mathbf{x}_S(n) - \mathbf{x}_A\right)
  39
  40 and the optimal position for the shells :math:`x_S(n+1)` thus follows
  41 from
  42
  43 .. 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
  44
  45 if we write
  46
  47 .. math:: \Delta \mathbf{x}_S = \mathbf{x}_S(n+1) - \mathbf{x}_S(n)
  48
  49 we finally obtain
  50
  51 .. math:: \Delta \mathbf{x}_S = \mathbf{F}_S/k_b
  52
  53 which then yields the algorithm to compute the next trial in the
  54 optimization of shell positions
  55
  56 .. math:: \mathbf{x}_S(n+1) ~=~ \mathbf{x}_S(n) + \mathbf{F}_S/k_b.