docs/reference-manual/algorithms/stochastic-dynamics.rst

   1 Stochastic Dynamics
   2 -------------------
   3
   4 Stochastic or velocity Langevin dynamics adds a friction and a noise
   5 term to Newton’s equations of motion, as
   6
   7 .. math::  m_i {{\mbox{d}}^2 \mathbf{r}_i \over {\mbox{d}}t^2} =
   8    - m_i \gamma_i {{\mbox{d}}\mathbf{r}_i \over {\mbox{d}}t} + \mathbf{F}_i(\mathbf{r}) + {\stackrel{\circ}{\mathbf{r}}}_i,
   9    :label: eqnSDeq
  10
  11 where :math:`\gamma_i` is the friction constant :math:`[1/\mbox{ps}]`
  12 and :math:`{\stackrel{\circ}{\mathbf{r}}}_i\!\!(t)` is a
  13 noise process with
  14 :math:`\langle {\stackrel{\circ}{r}}_i\!\!(t) {\stackrel{\circ}{r}}_j\!\!(t+s) \rangle = 2 m_i \gamma_i k_B T \delta(s) \delta_{ij}`. When :math:`1/\gamma_i`
  15 is large compared to the time scales present in the system, one could
  16 see stochastic dynamics as molecular dynamics with stochastic
  17 temperature-coupling. But any processes that take longer than
  18 :math:`1/\gamma_i`, e.g. hydrodynamics, will be dampened. Since each
  19 degree of freedom is coupled independently to a heat bath, equilibration
  20 of fast modes occurs rapidly. For simulating a system in vacuum there is
  21 the additional advantage that there is no accumulation of errors for the
  22 overall translational and rotational degrees of freedom. When
  23 :math:`1/\gamma_i` is small compared to the time scales present in the
  24 system, the dynamics will be completely different from MD, but the
  25 sampling is still correct.
  26
  27 In |Gromacs| there is one simple and efficient implementation. Its
  28 accuracy is equivalent to the normal MD leap-frog and Velocity Verlet
  29 integrator. It is nearly identical to the common way of discretizing the
  30 Langevin equation, but the friction and velocity term are applied in an
  31 impulse fashion \ :ref:`51 <refGoga2012>`. It can be described as:
  32
  33 .. math::  \begin{aligned}
  34    \mathbf{v}'  &~=~&   \mathbf{v}(t-{{\frac{1}{2}}{{\Delta t}}}) + \frac{1}{m}\mathbf{F}(t){{\Delta t}}\\
  35    \Delta\mathbf{v}     &~=~&   -\alpha \, \mathbf{v}'(t+{{\frac{1}{2}}{{\Delta t}}}) + \sqrt{\frac{k_B T}{m}(1 - \alpha^2)} \, {\mathbf{r}^G}_i \\
  36    \mathbf{r}(t+{{\Delta t}})   &~=~&   \mathbf{r}(t)+\left(\mathbf{v}' +\frac{1}{2}\Delta \mathbf{v}\right){{\Delta t}}
  37   \end{aligned}
  38   :label: eqnsd1int
  39
  40 .. math:: \begin{aligned}
  41    \mathbf{v}(t+{{\frac{1}{2}}{{\Delta t}}})  &~=~&   \mathbf{v}' + \Delta \mathbf{v} \\
  42    \alpha &~=~& 1 - e^{-\gamma {{\Delta t}}}\end{aligned}
  43    :label: eqnsd1xupd
  44
  45 where :math:`{\mathbf{r}^G}_i` is Gaussian distributed
  46 noise with :math:`\mu = 0`, :math:`\sigma = 1`. The velocity is first
  47 updated a full time step without friction and noise to get
  48 :math:`\mathbf{v}'`, identical to the normal update in
  49 leap-frog. The friction and noise are then applied as an impulse at step
  50 :math:`t+{{\Delta t}}`. The advantage of this scheme is that the
  51 velocity-dependent terms act at the full time step, which makes the
  52 correct integration of forces that depend on both coordinates and
  53 velocities, such as constraints and dissipative particle dynamics (DPD,
  54 not implented yet), straightforward. With constraints, the coordinate
  55 update :eq:`eqn. %s <eqnsd1xupd>` is split into a normal leap-frog update
  56 and a :math:`\Delta \mathbf{v}`. After both of these
  57 updates the constraints are applied to coordinates and velocities.
  58
  59 When using SD as a thermostat, an appropriate value for :math:`\gamma`
  60 is e.g. 0.5 ps\ :math:`^{-1}`, since this results in a friction that is
  61 lower than the internal friction of water, while it still provides
  62 efficient thermostatting.