:math:`T`:
.. math:: p(v_i) = \sqrt{\frac{m_i}{2 \pi kT}}\exp\left(-\frac{m_i v_i^2}{2kT}\right)
+ :label: eqnmaxwellboltzman
where :math:`k` is Boltzmann’s constant (see chapter :ref:`defunits`). To
accomplish this, normally distributed random numbers are generated by
density :math:`\rho_2`, when the inter-particle distance changes from
:math:`r_0` to :math:`r_t`, as:
-.. math::
-
- \langle \Delta V \rangle =
- \int_{0}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 V(r_t) G\!\left(\frac{r_t-r_0}{\sigma}\right) d r_0\, d r_t
+.. math:: \langle \Delta V \rangle =
+ \int_{0}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 V(r_t) G\!\left(\frac{r_t-r_0}{\sigma}\right) d r_0\, d r_t
+ :label: eqnverletbufenergy
To evaluate this analytically, we need to make some approximations.
First we replace :math:`V(r_t)` by a Taylor expansion around
:math:`r_c`, then we can move the lower bound of the integral over
:math:`r_0` to :math:`-\infty` which will simplify the result:
-.. math::
-
- \begin{aligned}
- \langle \Delta V \rangle &\approx&
- \int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[ V'(r_c) (r_t - r_c) +
- \nonumber\\
- & &
- \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
- V''(r_c)\frac{1}{2}(r_t - r_c)^2 +
- \nonumber\\
- & &
- \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
- V'''(r_c)\frac{1}{6}(r_t - r_c)^3 +
- \nonumber\\
- & &
- \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
- O \! \left((r_t - r_c)^4 \right)\Big] G\!\left(\frac{r_t-r_0}{\sigma}\right) d r_0 \, d r_t\end{aligned}
+.. math:: \begin{aligned}
+ \langle \Delta V \rangle &\approx&
+ \int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[ V'(r_c) (r_t - r_c) +
+ \nonumber\\
+ & &
+ \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
+ V''(r_c)\frac{1}{2}(r_t - r_c)^2 +
+ \nonumber\\
+ & &
+ \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
+ V'''(r_c)\frac{1}{6}(r_t - r_c)^3 +
+ \nonumber\\
+ & &
+ \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
+ O \! \left((r_t - r_c)^4 \right)\Big] G\!\left(\frac{r_t-r_0}{\sigma}\right) d r_0 \, d r_t\end{aligned}
+ :label: eqnverletaylor
Replacing the factor :math:`r_0^2` by :math:`(r_\ell + \sigma)^2`,
which results in a slight overestimate, allows us to calculate the
integrals analytically:
-.. math::
-
- \begin{aligned}
- \langle \Delta V \rangle \!
- &\approx&
- 4 \pi (r_\ell+\sigma)^2 \rho_2
- \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[ V'(r_c) (r_t - r_c) +
- \nonumber\\
- & &
- \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
- V''(r_c)\frac{1}{2}(r_t - r_c)^2 +
- \nonumber\\
- & &
- \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
- V'''(r_c)\frac{1}{6}(r_t - r_c)^3 \Big] G\!\left(\frac{r_t-r_0}{\sigma}\right)
- d r_0 \, d r_t\\
- &=&
- 4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{
- \frac{1}{2}V'(r_c)\left[r_b \sigma G\!\left(\frac{r_b}{\sigma}\right) - (r_b^2+\sigma^2)E\!\left(\frac{r_b}{\sigma}\right) \right] +
- \nonumber\\
- & &
- \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
- \frac{1}{6}V''(r_c)\left[ \sigma(r_b^2+2\sigma^2) G\!\left(\frac{r_b}{\sigma}\right) - r_b(r_b^2+3\sigma^2 ) E\!\left(\frac{r_b}{\sigma}\right) \right] +
- \nonumber\\
- & &
- \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
- \frac{1}{24}V'''(r_c)\bigg[ r_b\sigma(r_b^2+5\sigma^2) G\!\left(\frac{r_b}{\sigma}\right)
- \nonumber\\
- & &
- \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ \frac{1}{24}V'''(r_c)\bigg[ }
- - (r_b^4+6r_b^2\sigma^2+3\sigma^4 ) E\!\left(\frac{r_b}{\sigma}\right) \bigg]
- \bigg\}\end{aligned}
+.. math:: \begin{aligned}
+ \langle \Delta V \rangle \!
+ &\approx&
+ 4 \pi (r_\ell+\sigma)^2 \rho_2
+ \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[ V'(r_c) (r_t - r_c) +
+ \nonumber\\
+ & &
+ \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
+ V''(r_c)\frac{1}{2}(r_t - r_c)^2 +
+ \nonumber\\
+ & &
+ \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
+ V'''(r_c)\frac{1}{6}(r_t - r_c)^3 \Big] G\!\left(\frac{r_t-r_0}{\sigma}\right)
+ d r_0 \, d r_t\\
+ &=&
+ 4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{
+ \frac{1}{2}V'(r_c)\left[r_b \sigma G\!\left(\frac{r_b}{\sigma}\right) - (r_b^2+\sigma^2)E\!\left(\frac{r_b}{\sigma}\right) \right] +
+ \nonumber\\
+ & &
+ \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
+ \frac{1}{6}V''(r_c)\left[ \sigma(r_b^2+2\sigma^2) G\!\left(\frac{r_b}{\sigma}\right) - r_b(r_b^2+3\sigma^2 ) E\!\left(\frac{r_b}{\sigma}\right) \right] +
+ \nonumber\\
+ & &
+ \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
+ \frac{1}{24}V'''(r_c)\bigg[ r_b\sigma(r_b^2+5\sigma^2) G\!\left(\frac{r_b}{\sigma}\right)
+ \nonumber\\
+ & &
+ \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ \frac{1}{24}V'''(r_c)\bigg[ }
+ - (r_b^4+6r_b^2\sigma^2+3\sigma^4 ) E\!\left(\frac{r_b}{\sigma}\right) \bigg]
+ \bigg\}\end{aligned}
+ :label: eqnverletanalytical
where :math:`G(x)` is a Gaussian distribution with 0 mean and unit
variance and :math:`E(x)=\frac{1}{2}\mathrm{erfc}(x/\sqrt{2})`. We
:math:`{\mathbf{r}_{ij}}` connecting images within the
cut-off :math:`R_c` can be found by constructing:
-.. math::
-
- \begin{aligned}
- \mathbf{r}''' & = & \mathbf{r}_j-\mathbf{r}_i \\
- \mathbf{r}'' & = & \mathbf{r}''' - \mathbf{c}*\mathrm{round}(r'''_z/c_z) \\
- \mathbf{r}' & = & \mathbf{r}'' - \mathbf{b}*\mathrm{round}(r''_y/b_y) \\
- \mathbf{r}_{ij} & = & \mathbf{r}' - \mathbf{a}*\mathrm{round}(r'_x/a_x)
- \end{aligned}
+.. math:: \begin{aligned}
+ \mathbf{r}''' & = & \mathbf{r}_j-\mathbf{r}_i \\
+ \mathbf{r}'' & = & \mathbf{r}''' - \mathbf{c}*\mathrm{round}(r'''_z/c_z) \\
+ \mathbf{r}' & = & \mathbf{r}'' - \mathbf{b}*\mathrm{round}(r''_y/b_y) \\
+ \mathbf{r}_{ij} & = & \mathbf{r}' - \mathbf{a}*\mathrm{round}(r'_x/a_x)
+ \end{aligned}
+ :label: eqnsearchvec
When distances between two particles in a triclinic box are needed that
do not obey :eq:`eqn. %s <eqnboxrot>`, many shifts of
:math:`N`-particle system:
.. math:: E_{kin} = {\frac{1}{2}}\sum_{i=1}^N m_i v_i^2
+ :label: eqntempEkin
From this the absolute temperature :math:`T` can be computed using:
number of degrees of freedom which can be computed from:
.. math:: N_{\mathrm{df}} ~=~ 3 N - N_c - N_{\mathrm{com}}
+ :label: eqndofcoupling
Here :math:`N_c` is the number of *constraints* imposed on the system.
When performing molecular dynamics :math:`N_{\mathrm{com}}=3` additional
for group :math:`i` is:
.. math:: N^i_{\mathrm{df}} ~=~ (3 N^i - N^i_c) \frac{3 N - N_c - N_{\mathrm{com}}}{3 N - N_c}
+ :label: eqndofonecouplgroup
The kinetic energy can also be written as a tensor, which is necessary
for pressure calculation in a triclinic system, or systems where shear
forces are imposed:
.. math:: {\bf E}_{\mathrm{kin}} = {\frac{1}{2}}\sum_i^N m_i {\mathbf{v}_i}\otimes {\mathbf{v}_i}
+ :label: eqnEkintensor
Pressure and virial
^^^^^^^^^^^^^^^^^^^
algorithm, whose position-update relation is
.. math:: \mathbf{r}(t+{{\Delta t}})~=~2\mathbf{r}(t) - \mathbf{r}(t-{{\Delta t}}) + \frac{1}{m}\mathbf{F}(t){{\Delta t}}^2+O({{\Delta t}}^4)
+ :label: eqnleapfrogp
The algorithm is of third order in :math:`\mathbf{r}` and
is time-reversible. See ref. \ :ref:`24 <refBerendsen86b>` for the
from time :math:`t = 0` to time :math:`t` by applying the evolution
operator
-.. math::
-
- \begin{aligned}
- \Gamma(t) &=& \exp(iLt) \Gamma(0) \nonumber \\
- iL &=& \dot{\Gamma}\cdot \nabla_{\Gamma},\end{aligned}
+.. math:: \begin{aligned}
+ \Gamma(t) &=& \exp(iLt) \Gamma(0) \nonumber \\
+ iL &=& \dot{\Gamma}\cdot \nabla_{\Gamma},\end{aligned}
+ :label: eqnevoluoperator
where :math:`L` is the Liouville operator, and :math:`\Gamma` is the
multidimensional vector of independent variables (positions and
succession to evolve the system as
.. math:: \Gamma(t) = \prod_{i=1}^P \exp(iL{{\Delta t}}) \Gamma(0)
+ :label: eqnevolvesystem
For NVE dynamics, the Liouville operator is
-.. math::
-
- \begin{aligned}
- iL = \sum_{i=1}^{N} {{\mathbf{v}}}_i \cdot \nabla_{{{\mathbf{r}}}_i} + \sum_{i=1}^N \frac{1}{m_i}{{\mathbf{F}}}(r_i) \cdot \nabla_{{{\mathbf{v}}}_i}.\end{aligned}
+.. math:: \begin{aligned}
+ iL = \sum_{i=1}^{N} {{\mathbf{v}}}_i \cdot \nabla_{{{\mathbf{r}}}_i} + \sum_{i=1}^N \frac{1}{m_i}{{\mathbf{F}}}(r_i) \cdot \nabla_{{{\mathbf{v}}}_i}.\end{aligned}
+ :label: eqnliouvilleoperator
This can be split into two additive operators
-.. math::
-
- \begin{aligned}
- iL_1 &=& \sum_{i=1}^N \frac{1}{m_i}{{\mathbf{F}}}(r_i) \cdot \nabla_{{{\mathbf{v}}}_i} \nonumber \\
- iL_2 &=& \sum_{i=1}^{N} {{\mathbf{v}}}_i \cdot \nabla_{{{\mathbf{r}}}_i} \end{aligned}
+.. math:: \begin{aligned}
+ iL_1 &=& \sum_{i=1}^N \frac{1}{m_i}{{\mathbf{F}}}(r_i) \cdot \nabla_{{{\mathbf{v}}}_i} \nonumber \\
+ iL_2 &=& \sum_{i=1}^{N} {{\mathbf{v}}}_i \cdot \nabla_{{{\mathbf{r}}}_i} \end{aligned}
+ :label: eqnlotwoadditive
Then a short-time, symmetric, and thus reversible approximation of the
true dynamics will be
:math:`{{\frac{1}{2}}{{\Delta t}}}` is the final velocity half step. For
future times :math:`t = n{{\Delta t}}`, this becomes
-.. math::
-
- \begin{aligned}
- \exp(iLn{{\Delta t}}) &\approx& \left(\exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}})\right)^n \nonumber \\
- &\approx& \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \bigg(\exp(iL_1{{\Delta t}}) \exp(iL_2{{\Delta t}})\bigg)^{n-1} \nonumber \\
- & & \;\;\;\; \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \end{aligned}
+.. math:: \begin{aligned}
+ \exp(iLn{{\Delta t}}) &\approx& \left(\exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}})\right)^n \nonumber \\
+ &\approx& \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \bigg(\exp(iL_1{{\Delta t}}) \exp(iL_2{{\Delta t}})\bigg)^{n-1} \nonumber \\
+ & & \;\;\;\; \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \end{aligned}
+ :label: eqntrottertimestep
This formalism allows us to easily see the difference between the
different flavors of Verlet integrators. The leap-frog integrator can be
:math:`\exp\left(iL_1 {\Delta t}\right)` term, instead of the half-step
velocity term, yielding
-.. math::
-
- \begin{aligned}
- \exp(iLn{\Delta t}) &=& \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {{\Delta t}}\right) + \mathcal{O}({{\Delta t}}^3).\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iLn{\Delta t}) &=& \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {{\Delta t}}\right) + \mathcal{O}({{\Delta t}}^3).\end{aligned}
+ :label: eqnleapfroghalfvel
Here, the full step in velocity is between
:math:`t-{{\frac{1}{2}}{{\Delta t}}}` and
velocity half steps in velocity Verlet. For future times
:math:`t = n{{\Delta t}}`, this becomes
-.. math::
-
- \begin{aligned}
- \exp(iLn{\Delta t}) &\approx& \bigg(\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {{\Delta t}}\right) \bigg)^{n}.\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iLn{\Delta t}) &\approx& \bigg(\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {{\Delta t}}\right) \bigg)^{n}.\end{aligned}
+ :label: eqnvelverlethalfvel
Although at first this does not appear symmetric, as long as the full
velocity step is between :math:`t-{{\frac{1}{2}}{{\Delta t}}}` and
thus the temperature only at time :math:`t`; the kinetic energy is then
a sum over all particles
-.. math::
-
- \begin{aligned}
- KE_{\mathrm{full}}(t) &=& \sum_i \left(\frac{1}{2m_i}\mathbf{v}_i(t)\right)^2 \nonumber\\
- &=& \sum_i \frac{1}{2m_i}\left(\frac{1}{2}\mathbf{v}_i(t-{{\frac{1}{2}}{{\Delta t}}})+\frac{1}{2}\mathbf{v}_i(t+{{\frac{1}{2}}{{\Delta t}}})\right)^2,\end{aligned}
+.. math:: \begin{aligned}
+ KE_{\mathrm{full}}(t) &=& \sum_i \left(\frac{1}{2m_i}\mathbf{v}_i(t)\right)^2 \nonumber\\
+ &=& \sum_i \frac{1}{2m_i}\left(\frac{1}{2}\mathbf{v}_i(t-{{\frac{1}{2}}{{\Delta t}}})+\frac{1}{2}\mathbf{v}_i(t+{{\frac{1}{2}}{{\Delta t}}})\right)^2,\end{aligned}
+ :label: eqnTrotterEkin
with the square on the *outside* of the average. Standard leap-frog
calculates the kinetic energy at time :math:`t` based on the average
kinetic energies at the timesteps :math:`t+{{\frac{1}{2}}{{\Delta t}}}`
and :math:`t-{{\frac{1}{2}}{{\Delta t}}}`, or the sum over all particles
-.. math::
-
- \begin{aligned}
- KE_{\mathrm{average}}(t) &=& \sum_i \frac{1}{2m_i}\left(\frac{1}{2}\mathbf{v}_i(t-{{\frac{1}{2}}{{\Delta t}}})^2+\frac{1}{2}\mathbf{v}_i(t+{{\frac{1}{2}}{{\Delta t}}})^2\right),\end{aligned}
+.. math:: \begin{aligned}
+ KE_{\mathrm{average}}(t) &=& \sum_i \frac{1}{2m_i}\left(\frac{1}{2}\mathbf{v}_i(t-{{\frac{1}{2}}{{\Delta t}}})^2+\frac{1}{2}\mathbf{v}_i(t+{{\frac{1}{2}}{{\Delta t}}})^2\right),\end{aligned}
+ :label: eqnTrottersumparticles
where the square is *inside* the average.
(:eq:`eqn. %s <eqnTcoupling>`):
.. math:: \tau = 2 C_V \tau_T / N_{df} k
+ :label: eqnTcoupltau
where :math:`C_V` is the total heat capacity of the system, :math:`k`
is Boltzmann’s constant, and :math:`N_{df}` is the total number of
The thermostat modifies the kinetic energy at each scaling step by:
.. math:: \Delta E_k = (\lambda - 1)^2 E_k
+ :label: eqnThermostat
The sum of these changes over the run needs to subtracted from the
total energy to obtain the conserved energy quantity.
where the equation of motion for the heat bath parameter :math:`\xi` is:
.. math:: \frac {{\mbox{d}}p_{\xi}}{{\mbox{d}}t} = \left( T - T_0 \right).
+ :label: eqnNHheatbath
The reference temperature is denoted :math:`T_0`, while :math:`T` is
the current instantaneous temperature of the system. The strength of the
The conserved quantity for the Nosé-Hoover equations of motion is not
the total energy, but rather
-.. math::
-
- \begin{aligned}
- H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) +\frac{p_{\xi}^2}{2Q} + N_fkT\xi,\end{aligned}
+.. math:: \begin{aligned}
+ H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) +\frac{p_{\xi}^2}{2Q} + N_fkT\xi,\end{aligned}
+ :label: eqnNHconservedbasic
where :math:`N_f` is the total number of degrees of freedom.
:math:`T_0` via:
.. math:: Q = \frac {\tau_T^2 T_0}{4 \pi^2}.
+ :label: eqnNHQ
This provides a much more intuitive way of selecting the Nosé-Hoover
coupling strength (similar to the weak-coupling relaxation), and in
The conserved quantity for Nosé-Hoover chains is
-.. math::
-
- \begin{aligned}
- H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) +\sum_{k=1}^M\frac{p^2_{\xi_k}}{2Q^{\prime}_k} + N_fkT\xi_1 + kT\sum_{k=2}^M \xi_k \end{aligned}
+.. math:: \begin{aligned}
+ H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) +\sum_{k=1}^M\frac{p^2_{\xi_k}}{2Q^{\prime}_k} + N_fkT\xi_1 + kT\sum_{k=2}^M \xi_k \end{aligned}
+ :label: eqnNHconservedquantity
The values and velocities of the Nosé-Hoover thermostat variables are
generally not included in the output, as they take up a fair amount of
operator as
.. math:: iL = iL_1 + iL_2 + iL_{\mathrm{NHC}},
+ :label: eqnNHTrotter
where
-.. math::
-
- \begin{aligned}
- iL_1 &=& \sum_{i=1}^N \left[\frac{{{\mathbf{p}}}_i}{m_i}\right]\cdot \frac{\partial}{\partial {{\mathbf{r}}}_i} \nonumber \\
- iL_2 &=& \sum_{i=1}^N {{\mathbf{F}}}_i\cdot \frac{\partial}{\partial {{\mathbf{p}}}_i} \nonumber \\
- iL_{\mathrm{NHC}} &=& \sum_{i=1}^N-\frac{p_{\xi}}{Q}{{\mathbf{v}}}_i\cdot \nabla_{{{\mathbf{v}}}_i} +\frac{p_{\xi}}{Q}\frac{\partial }{\partial \xi} + \left( T - T_0 \right)\frac{\partial }{\partial p_{\xi}}\end{aligned}
+.. math:: \begin{aligned}
+ iL_1 &=& \sum_{i=1}^N \left[\frac{{{\mathbf{p}}}_i}{m_i}\right]\cdot \frac{\partial}{\partial {{\mathbf{r}}}_i} \nonumber \\
+ iL_2 &=& \sum_{i=1}^N {{\mathbf{F}}}_i\cdot \frac{\partial}{\partial {{\mathbf{p}}}_i} \nonumber \\
+ iL_{\mathrm{NHC}} &=& \sum_{i=1}^N-\frac{p_{\xi}}{Q}{{\mathbf{v}}}_i\cdot \nabla_{{{\mathbf{v}}}_i} +\frac{p_{\xi}}{Q}\frac{\partial }{\partial \xi} + \left( T - T_0 \right)\frac{\partial }{\partial p_{\xi}}\end{aligned}
+ :label: eqnNHTrotter2
For standard velocity Verlet with Nosé-Hoover temperature control, this
becomes
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
- &&\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) + \mathcal{O}({{\Delta t}}^3).\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
+ &&\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) + \mathcal{O}({{\Delta t}}^3).\end{aligned}
+ :label: eqnNHTrotter3
For half-step-averaged temperature control using *md-vv-avek*, this
decomposition will not work, since we do not have the full step
:eq:`Eq. %s <eqhalfstepNHCintegrator>` just before the
:math:`\exp\left(iL_1 {\Delta t}\right)` term, yielding:
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &=& \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \\
- &&\exp\left(iL_2 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) + \mathcal{O}({{\Delta t}}^3)\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &=& \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \\
+ &&\exp\left(iL_2 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) + \mathcal{O}({{\Delta t}}^3)\end{aligned}
+ :label: eqnNHleapfrog
and then using some algebra tricks to solve for some quantities are
required before they are actually calculated \ :ref:`36 <refHolian95>`.
towards a given reference pressure :math:`{\bf P}_0` according to
.. math:: \frac{{\mbox{d}}{\bf P}}{{\mbox{d}}t} = \frac{{\bf P}_0-{\bf P}}{\tau_p}.
+ :label: eqnberendsenpressure
The scaling matrix :math:`\mu` is given by
.. math:: \mu_{ij}
- = \delta_{ij} - \frac{n_\mathrm{PC}\Delta t}{3\, \tau_p} \beta_{ij} \{P_{0ij} - P_{ij}(t) \}.
- :label: eqnmu
+ = \delta_{ij} - \frac{n_\mathrm{PC}\Delta t}{3\, \tau_p} \beta_{ij} \{P_{0ij} - P_{ij}(t) \}.
+ :label: eqnmu
Here, :math:`\beta` is the isothermal compressibility of the system. In
most cases this will be a diagonal matrix, with equal elements on the
scaling, which is usually less than :math:`10^{-4}`. The actual scaling
matrix :math:`\mu'` is
-.. math::
-
- \mathbf{\mu'} =
- \left(\begin{array}{ccc}
- \mu_{xx} & \mu_{xy} + \mu_{yx} & \mu_{xz} + \mu_{zx} \\
- 0 & \mu_{yy} & \mu_{yz} + \mu_{zy} \\
- 0 & 0 & \mu_{zz}
- \end{array}\right).
+.. math:: \mathbf{\mu'} =
+ \left(\begin{array}{ccc}
+ \mu_{xx} & \mu_{xy} + \mu_{yx} & \mu_{xz} + \mu_{zx} \\
+ 0 & \mu_{yy} & \mu_{yz} + \mu_{zy} \\
+ 0 & 0 & \mu_{zz}
+ \end{array}\right).
+ :label: eqnberendsenpressurescaling
The velocities are neither scaled nor rotated. Since the equations of
motion are modified by pressure coupling, the conserved energy quantity
the barostat applies to the system every step needs to be subtracted
from the total energy to obtain the conserved energy quantity:
-.. math::
-
- - \sum_{i,j} (\mu_{ij} -\delta_{ij}) P_{ij} V =
- \sum_{i,j} 2(\mu_{ij} -\delta_{ij}) \Xi_{ij}
+.. math:: - \sum_{i,j} (\mu_{ij} -\delta_{ij}) P_{ij} V =
+ \sum_{i,j} 2(\mu_{ij} -\delta_{ij}) \Xi_{ij}
+ :label: eqnberendsenpressureconserved
where :math:`\delta_{ij}` is the Kronecker delta and :math:`{\bf \Xi}`
is the virial. Note that the factor 2 originates from the factor
obey the matrix equation of motion [2]_
.. math:: \frac{{\mbox{d}}\mathbf{b}^2}{{\mbox{d}}t^2}= V \mathbf{W}^{-1} \mathbf{b}'^{-1} \left( \mathbf{P} - \mathbf{P}_{ref}\right).
+ :label: eqnPRpressure
The volume of the box is denoted :math:`V`, and
:math:`\mathbf{W}` is a matrix parameter that determines
it simple we only show the Parrinello-Rahman modification here. The
modified Hamiltonian, which will be conserved, is:
-.. math::
-
- E_\mathrm{pot} + E_\mathrm{kin} + \sum_i P_{ii} V +
- \sum_{i,j} \frac{1}{2} W_{ij} \left( \frac{{\mbox{d}}b_{ij}}{{\mbox{d}}t} \right)^2
+.. math:: E_\mathrm{pot} + E_\mathrm{kin} + \sum_i P_{ii} V +
+ \sum_{i,j} \frac{1}{2} W_{ij} \left( \frac{{\mbox{d}}b_{ij}}{{\mbox{d}}t} \right)^2
+ :label: eqnPRpressureconserved
The equations of motion for the atoms, obtained from the Hamiltonian
are:
-.. math::
-
- \begin{aligned}
- \frac {{\mbox{d}}^2\mathbf{r}_i}{{\mbox{d}}t^2} & = & \frac{\mathbf{F}_i}{m_i} -
- \mathbf{M} \frac{{\mbox{d}}\mathbf{r}_i}{{\mbox{d}}t} , \\ \mathbf{M} & = & \mathbf{b}^{-1} \left[
- \mathbf{b} \frac{{\mbox{d}}\mathbf{b}'}{{\mbox{d}}t} + \frac{{\mbox{d}}\mathbf{b}}{{\mbox{d}}t} \mathbf{b}'
- \right] \mathbf{b}'^{-1}.
- \end{aligned}
+.. math:: \begin{aligned}
+ \frac {{\mbox{d}}^2\mathbf{r}_i}{{\mbox{d}}t^2} & = & \frac{\mathbf{F}_i}{m_i} -
+ \mathbf{M} \frac{{\mbox{d}}\mathbf{r}_i}{{\mbox{d}}t} , \\ \mathbf{M} & = & \mathbf{b}^{-1} \left[
+ \mathbf{b} \frac{{\mbox{d}}\mathbf{b}'}{{\mbox{d}}t} + \frac{{\mbox{d}}\mathbf{b}}{{\mbox{d}}t} \mathbf{b}'
+ \right] \mathbf{b}'^{-1}.
+ \end{aligned}
+ :label: eqnPRpressuremotion
This extra term has the appearance of a friction, but it should be
noted that it is ficticious, and rather an effect of the
time constant :math:`\tau_p` in the input file (:math:`L` is the largest
box matrix element):
-.. math::
-
- \left(
- \mathbf{W}^{-1} \right)_{ij} = \frac{4 \pi^2 \beta_{ij}}{3 \tau_p^2 L}.
+.. math:: \left(
+ \mathbf{W}^{-1} \right)_{ij} = \frac{4 \pi^2 \beta_{ij}}{3 \tau_p^2 L}.
+ :label: eqnPRpressuretimeconst
Just as for the Nosé-Hoover thermostat, you should realize that the
Parrinello-Rahman time constant is *not* equivalent to the relaxation
be calculated from the difference between the normal and the lateral
pressure
-.. math::
-
- \begin{aligned}
- \gamma(t) & = &
- \frac{1}{n} \int_0^{L_z}
- \left\{ P_{zz}(z,t) - \frac{P_{xx}(z,t) + P_{yy}(z,t)}{2} \right\} \mbox{d}z \\
- & = &
- \frac{L_z}{n} \left\{ P_{zz}(t) - \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\},\end{aligned}
+.. math:: \begin{aligned}
+ \gamma(t) & = &
+ \frac{1}{n} \int_0^{L_z}
+ \left\{ P_{zz}(z,t) - \frac{P_{xx}(z,t) + P_{yy}(z,t)}{2} \right\} \mbox{d}z \\
+ & = &
+ \frac{L_z}{n} \left\{ P_{zz}(t) - \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\},\end{aligned}
+ :label: eqnsurftenscoupl
where :math:`L_z` is the height of the box and :math:`n` is the number
of surfaces. The pressure in the z-direction is corrected by scaling the
height of the box with :math:`\mu_{zz}`
.. math:: \Delta P_{zz} = \frac{\Delta t}{\tau_p} \{ P_{0zz} - P_{zz}(t) \}
+ :label: eqnzpressure
.. math:: \mu_{zz} = 1 + \beta_{zz} \Delta P_{zz}
+ :label: eqnzpressure2
This is similar to normal pressure coupling, except that the factor of
:math:`1/3` is missing. The pressure correction in the
surface tension to the reference value :math:`\gamma_0`. The correction
factor for the box length in the :math:`x`/:math:`y`-direction is
-.. math::
-
- \mu_{x/y} = 1 + \frac{\Delta t}{2\,\tau_p} \beta_{x/y}
- \left( \frac{n \gamma_0}{\mu_{zz} L_z}
- - \left\{ P_{zz}(t)+\Delta P_{zz} - \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\}
- \right)
+.. math:: \mu_{x/y} = 1 + \frac{\Delta t}{2\,\tau_p} \beta_{x/y}
+ \left( \frac{n \gamma_0}{\mu_{zz} L_z}
+ - \left\{ P_{zz}(t)+\Delta P_{zz} - \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\}
+ \right)
+ :label: eqnboxlengthcorr
The value of :math:`\beta_{zz}` is more critical than with normal
pressure coupling. Normally an incorrect compressibility will just scale
The isobaric equations are
-.. math::
-
- \begin{aligned}
- \dot{{{\mathbf{r}}}}_i &=& \frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W} {{\mathbf{r}}}_i \nonumber \\
- \frac{\dot{{{\mathbf{p}}}}_i}{m_i} &=& \frac{1}{m_i}{{\mathbf{F}}}_i - \alpha\frac{{p_{\epsilon}}}{W} \frac{{{\mathbf{p}}}_i}{m_i} \nonumber \\
- \dot{\epsilon} &=& \frac{{p_{\epsilon}}}{W} \nonumber \\
- \frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}(P_{\mathrm{int}} - P) + (\alpha-1)\left(\sum_{n=1}^N\frac{{{\mathbf{p}}}_i^2}{m_i}\right),\\\end{aligned}
+.. math:: \begin{aligned}
+ \dot{{{\mathbf{r}}}}_i &=& \frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W} {{\mathbf{r}}}_i \nonumber \\
+ \frac{\dot{{{\mathbf{p}}}}_i}{m_i} &=& \frac{1}{m_i}{{\mathbf{F}}}_i - \alpha\frac{{p_{\epsilon}}}{W} \frac{{{\mathbf{p}}}_i}{m_i} \nonumber \\
+ \dot{\epsilon} &=& \frac{{p_{\epsilon}}}{W} \nonumber \\
+ \frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}(P_{\mathrm{int}} - P) + (\alpha-1)\left(\sum_{n=1}^N\frac{{{\mathbf{p}}}_i^2}{m_i}\right),\\\end{aligned}
+ :label: eqnMTTKisobaric
where
-.. math::
-
- \begin{aligned}
- P_{\mathrm{int}} &=& P_{\mathrm{kin}} -P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{{{\mathbf{p}}}_i^2}{2m_i} - {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\
- \right)\right].\end{aligned}
+.. math:: \begin{aligned}
+ P_{\mathrm{int}} &=& P_{\mathrm{kin}} -P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{{{\mathbf{p}}}_i^2}{2m_i} - {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\
+ \right)\right].\end{aligned}
+ :label: eqnMTTKisobaric2
The terms including :math:`\alpha` are required to make phase space
incompressible \ :ref:`41 <refTuckerman2006>`. The :math:`\epsilon`
acceleration term can be rewritten as
-.. math::
-
- \begin{aligned}
- \frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}\left(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P\right)\end{aligned}
+.. math:: \begin{aligned}
+ \frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}\left(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P\right)\end{aligned}
+ :label: eqnMTTKaccel
In terms of velocities, these equations become
-.. math::
-
- \begin{aligned}
- \dot{{{\mathbf{r}}}}_i &=& {{\mathbf{v}}}_i + {v_{\epsilon}}{{\mathbf{r}}}_i \nonumber \\
- \dot{{{\mathbf{v}}}}_i &=& \frac{1}{m_i}{{\mathbf{F}}}_i - \alpha{v_{\epsilon}}{{\mathbf{v}}}_i \nonumber \\
- \dot{\epsilon} &=& {v_{\epsilon}}\nonumber \\
- \dot{{v_{\epsilon}}} &=& \frac{3V}{W}(P_{\mathrm{int}} - P) + (\alpha-1)\left( \sum_{n=1}^N \frac{1}{2} m_i {{\mathbf{v}}}_i^2\right)\nonumber \\
- P_{\mathrm{int}} &=& P_{\mathrm{kin}} - P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{1}{2} m_i{{\mathbf{v}}}_i^2 - {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\right)\right]\end{aligned}
+.. math:: \begin{aligned}
+ \dot{{{\mathbf{r}}}}_i &=& {{\mathbf{v}}}_i + {v_{\epsilon}}{{\mathbf{r}}}_i \nonumber \\
+ \dot{{{\mathbf{v}}}}_i &=& \frac{1}{m_i}{{\mathbf{F}}}_i - \alpha{v_{\epsilon}}{{\mathbf{v}}}_i \nonumber \\
+ \dot{\epsilon} &=& {v_{\epsilon}}\nonumber \\
+ \dot{{v_{\epsilon}}} &=& \frac{3V}{W}(P_{\mathrm{int}} - P) + (\alpha-1)\left( \sum_{n=1}^N \frac{1}{2} m_i {{\mathbf{v}}}_i^2\right)\nonumber \\
+ P_{\mathrm{int}} &=& P_{\mathrm{kin}} - P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{1}{2} m_i{{\mathbf{v}}}_i^2 - {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\right)\right]\end{aligned}
+ :label: eqnMTTKvel
For these equations, the conserved quantity is
-.. math::
-
- \begin{aligned}
- H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i^2}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) + \frac{p_\epsilon}{2W} + PV\end{aligned}
+.. math:: \begin{aligned}
+ H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i^2}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) + \frac{p_\epsilon}{2W} + PV\end{aligned}
+ :label: eqnMTTKconserved
The next step is to add temperature control. Adding Nosé-Hoover chains,
including to the barostat degree of freedom, where we use :math:`\eta`
for the barostat Nosé-Hoover variables, and :math:`Q^{\prime}` for the
coupling constants of the thermostats of the barostats, we get
-.. math::
-
- \begin{aligned}
- \dot{{{\mathbf{r}}}}_i &=& \frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W} {{\mathbf{r}}}_i \nonumber \\
- \frac{\dot{{{\mathbf{p}}}}_i}{m_i} &=& \frac{1}{m_i}{{\mathbf{F}}}_i - \alpha\frac{{p_{\epsilon}}}{W} \frac{{{\mathbf{p}}}_i}{m_i} - \frac{p_{\xi_1}}{Q_1}\frac{{{\mathbf{p}}}_i}{m_i}\nonumber \\
- \dot{\epsilon} &=& \frac{{p_{\epsilon}}}{W} \nonumber \\
- \frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P) -\frac{p_{\eta_1}}{Q^{\prime}_1}{p_{\epsilon}}\nonumber \\
- \dot{\xi}_k &=& \frac{p_{\xi_k}}{Q_k} \nonumber \\
- \dot{\eta}_k &=& \frac{p_{\eta_k}}{Q^{\prime}_k} \nonumber \\
- \dot{p}_{\xi_k} &=& G_k - \frac{p_{\xi_{k+1}}}{Q_{k+1}} \;\;\;\; k=1,\ldots, M-1 \nonumber \\
- \dot{p}_{\eta_k} &=& G^\prime_k - \frac{p_{\eta_{k+1}}}{Q^\prime_{k+1}} \;\;\;\; k=1,\ldots, M-1 \nonumber \\
- \dot{p}_{\xi_M} &=& G_M \nonumber \\
- \dot{p}_{\eta_M} &=& G^\prime_M, \nonumber \\\end{aligned}
+.. math:: \begin{aligned}
+ \dot{{{\mathbf{r}}}}_i &=& \frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W} {{\mathbf{r}}}_i \nonumber \\
+ \frac{\dot{{{\mathbf{p}}}}_i}{m_i} &=& \frac{1}{m_i}{{\mathbf{F}}}_i - \alpha\frac{{p_{\epsilon}}}{W} \frac{{{\mathbf{p}}}_i}{m_i} - \frac{p_{\xi_1}}{Q_1}\frac{{{\mathbf{p}}}_i}{m_i}\nonumber \\
+ \dot{\epsilon} &=& \frac{{p_{\epsilon}}}{W} \nonumber \\
+ \frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P) -\frac{p_{\eta_1}}{Q^{\prime}_1}{p_{\epsilon}}\nonumber \\
+ \dot{\xi}_k &=& \frac{p_{\xi_k}}{Q_k} \nonumber \\
+ \dot{\eta}_k &=& \frac{p_{\eta_k}}{Q^{\prime}_k} \nonumber \\
+ \dot{p}_{\xi_k} &=& G_k - \frac{p_{\xi_{k+1}}}{Q_{k+1}} \;\;\;\; k=1,\ldots, M-1 \nonumber \\
+ \dot{p}_{\eta_k} &=& G^\prime_k - \frac{p_{\eta_{k+1}}}{Q^\prime_{k+1}} \;\;\;\; k=1,\ldots, M-1 \nonumber \\
+ \dot{p}_{\xi_M} &=& G_M \nonumber \\
+ \dot{p}_{\eta_M} &=& G^\prime_M, \nonumber \\\end{aligned}
+ :label: eqnMTTKthermandbar
where
-.. math::
-
- \begin{aligned}
- P_{\mathrm{int}} &=& P_{\mathrm{kin}} - P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{{{\mathbf{p}}}_i^2}{2m_i} - {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\right)\right] \nonumber \\
- G_1 &=& \sum_{i=1}^N \frac{{{\mathbf{p}}}^2_i}{m_i} - N_f kT \nonumber \\
- G_k &=& \frac{p^2_{\xi_{k-1}}}{2Q_{k-1}} - kT \;\; k = 2,\ldots,M \nonumber \\
- G^\prime_1 &=& \frac{{p_{\epsilon}}^2}{2W} - kT \nonumber \\
- G^\prime_k &=& \frac{p^2_{\eta_{k-1}}}{2Q^\prime_{k-1}} - kT \;\; k = 2,\ldots,M\end{aligned}
+.. math:: \begin{aligned}
+ P_{\mathrm{int}} &=& P_{\mathrm{kin}} - P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{{{\mathbf{p}}}_i^2}{2m_i} - {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\right)\right] \nonumber \\
+ G_1 &=& \sum_{i=1}^N \frac{{{\mathbf{p}}}^2_i}{m_i} - N_f kT \nonumber \\
+ G_k &=& \frac{p^2_{\xi_{k-1}}}{2Q_{k-1}} - kT \;\; k = 2,\ldots,M \nonumber \\
+ G^\prime_1 &=& \frac{{p_{\epsilon}}^2}{2W} - kT \nonumber \\
+ G^\prime_k &=& \frac{p^2_{\eta_{k-1}}}{2Q^\prime_{k-1}} - kT \;\; k = 2,\ldots,M\end{aligned}
+ :label: eqnMTTKthermandbar2
The conserved quantity is now
-.. math::
-
- \begin{aligned}
- H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) + \frac{p^2_\epsilon}{2W} + PV + \nonumber \\
- \sum_{k=1}^M\frac{p^2_{\xi_k}}{2Q_k} +\sum_{k=1}^M\frac{p^2_{\eta_k}}{2Q^{\prime}_k} + N_fkT\xi_1 + kT\sum_{i=2}^M \xi_k + kT\sum_{k=1}^M \eta_k\end{aligned}
+.. math:: \begin{aligned}
+ H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) + \frac{p^2_\epsilon}{2W} + PV + \nonumber \\
+ \sum_{k=1}^M\frac{p^2_{\xi_k}}{2Q_k} +\sum_{k=1}^M\frac{p^2_{\eta_k}}{2Q^{\prime}_k} + N_fkT\xi_1 + kT\sum_{i=2}^M \xi_k + kT\sum_{k=1}^M \eta_k\end{aligned}
+ :label: eqnMTTKthermandbarconserved
Returning to the Trotter decomposition formalism, for pressure control
and temperature control \ :ref:`35 <refMartyna1996>` we get:
-.. math::
-
- \begin{aligned}
- iL = iL_1 + iL_2 + iL_{\epsilon,1} + iL_{\epsilon,2} + iL_{\mathrm{NHC-baro}} + iL_{\mathrm{NHC}}\end{aligned}
+.. math:: \begin{aligned}
+ iL = iL_1 + iL_2 + iL_{\epsilon,1} + iL_{\epsilon,2} + iL_{\mathrm{NHC-baro}} + iL_{\mathrm{NHC}}\end{aligned}
+ :label: eqnMTTKthermandbarTrotter
where “NHC-baro” corresponds to the Nosè-Hoover chain of the barostat,
and NHC corresponds to the NHC of the particles,
-.. math::
-
- \begin{aligned}
- iL_1 &=& \sum_{i=1}^N \left[\frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W}{{\mathbf{r}}}_i\right]\cdot \frac{\partial}{\partial {{\mathbf{r}}}_i} \\
- iL_2 &=& \sum_{i=1}^N {{\mathbf{F}}}_i - \alpha \frac{{p_{\epsilon}}}{W}{{\mathbf{p}}}_i \cdot \frac{\partial}{\partial {{\mathbf{p}}}_i} \\
- iL_{\epsilon,1} &=& \frac{p_\epsilon}{W} \frac{\partial}{\partial \epsilon}\\
- iL_{\epsilon,2} &=& G_{\epsilon} \frac{\partial}{\partial p_\epsilon}\end{aligned}
+.. math:: \begin{aligned}
+ iL_1 &=& \sum_{i=1}^N \left[\frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W}{{\mathbf{r}}}_i\right]\cdot \frac{\partial}{\partial {{\mathbf{r}}}_i} \\
+ iL_2 &=& \sum_{i=1}^N {{\mathbf{F}}}_i - \alpha \frac{{p_{\epsilon}}}{W}{{\mathbf{p}}}_i \cdot \frac{\partial}{\partial {{\mathbf{p}}}_i} \\
+ iL_{\epsilon,1} &=& \frac{p_\epsilon}{W} \frac{\partial}{\partial \epsilon}\\
+ iL_{\epsilon,2} &=& G_{\epsilon} \frac{\partial}{\partial p_\epsilon}\end{aligned}
+ :label: eqnMTTKthermandbarTrotter2
and where
-.. math::
-
- \begin{aligned}
- G_{\epsilon} = 3V\left(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P\right)\end{aligned}
+.. math:: \begin{aligned}
+ G_{\epsilon} = 3V\left(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P\right)\end{aligned}
+ :label: eqnMTTKthermandbarTrotter3
Using the Trotter decomposition, we get
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right)\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \nonumber \\
- &&\exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \nonumber \\
- &&\exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \nonumber \nonumber \\
- &&\exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \nonumber \nonumber \\
- &&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right)\exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right) + \mathcal{O}({\Delta t}^3)\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right)\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \nonumber \\
+ &&\exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \nonumber \\
+ &&\exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \nonumber \nonumber \\
+ &&\exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \nonumber \nonumber \\
+ &&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right)\exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right) + \mathcal{O}({\Delta t}^3)\end{aligned}
+ :label: eqnMTTKthermandbarTrotterdecomp
The action of :math:`\exp\left(iL_1 {\Delta t}\right)` comes from the
solution of the the differential equation
:math:`t=\Delta t`. This yields the evolution
.. math:: {{\mathbf{r}}}_i({\Delta t}) = {{\mathbf{r}}}_i(0)e^{{v_{\epsilon}}{\Delta t}} + \Delta t {{\mathbf{v}}}_i(0) e^{{v_{\epsilon}}{\Delta t}/2} {\frac{\sinh{\left( {v_{\epsilon}}{\Delta t}/2\right)}}{{v_{\epsilon}}{\Delta t}/2}}.
+ :label: eqnMTTKthermandbarTrotterevol
The action of :math:`\exp\left(iL_2 {\Delta t}/2\right)` comes from the
solution of the differential equation
\alpha{v_{\epsilon}}{{\mathbf{v}}}_i`, yielding
.. math:: {{\mathbf{v}}}_i({\Delta t}/2) = {{\mathbf{v}}}_i(0)e^{-\alpha{v_{\epsilon}}{\Delta t}/2} + \frac{\Delta t}{2m_i}{{\mathbf{F}}}_i(0) e^{-\alpha{v_{\epsilon}}{\Delta t}/4}{\frac{\sinh{\left( \alpha{v_{\epsilon}}{\Delta t}/4\right)}}{\alpha{v_{\epsilon}}{\Delta t}/4}}.
+ :label: eqnMTTKthermandbarTrotterevol2
*md-vv-avek* uses the full step kinetic energies for determining the
pressure with the pressure control, but the half-step-averaged kinetic
energy for the temperatures, which can be written as a Trotter
decomposition as
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right)\nonumber \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
- &&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \\
- &&\exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right) + \mathcal{O}({\Delta t}^3)\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right)\nonumber \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
+ &&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \\
+ &&\exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC-baro}}{\Delta t}/2\right) + \mathcal{O}({\Delta t}^3)\end{aligned}
+ :label: eqnvvavekTrotterdecomp
With constraints, the equations become significantly more complicated,
in that each of these equations need to be solved iteratively for the
Standard velocity Verlet with Nosé-Hoover temperature control has a
Trotter expansion
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &\approx& \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
- &&\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right).\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &\approx& \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
+ &&\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right).\end{aligned}
+ :label:eqnVVNHTrotter
If the Nosé-Hoover chain is sufficiently slow with respect to the
motions of the system, we can write an alternate integrator over
:math:`n` steps for velocity Verlet as
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &\approx& (\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right)\left[\exp\left(iL_2 {\Delta t}/2\right)\right. \nonumber \\
- &&\left.\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right)\right]^n \exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right).\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &\approx& (\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right)\left[\exp\left(iL_2 {\Delta t}/2\right)\right. \nonumber \\
+ &&\left.\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right)\right]^n \exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right).\end{aligned}
+ :label: eqnVVNHTrotter2
For pressure control, this becomes
-.. math::
-
- \begin{aligned}
- \exp(iL{\Delta t}) &\approx& \exp\left(iL_{\mathrm{NHC-baro}}(n{\Delta t}/2)\right)\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right) \nonumber \nonumber \\
- &&\exp\left(iL_{\epsilon,2}(n{\Delta t}/2)\right) \left[\exp\left(iL_2 {\Delta t}/2\right)\right. \nonumber \nonumber \\
- &&\exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \nonumber \nonumber \\
- &&\left.\exp\left(iL_2 {\Delta t}/2\right)\right]^n \exp\left(iL_{\epsilon,2}(n{\Delta t}/2)\right) \nonumber \nonumber \\
- &&\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right)\exp\left(iL_{\mathrm{NHC-baro}}(n{\Delta t}/2)\right),\end{aligned}
+.. math:: \begin{aligned}
+ \exp(iL{\Delta t}) &\approx& \exp\left(iL_{\mathrm{NHC-baro}}(n{\Delta t}/2)\right)\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right) \nonumber \nonumber \\
+ &&\exp\left(iL_{\epsilon,2}(n{\Delta t}/2)\right) \left[\exp\left(iL_2 {\Delta t}/2\right)\right. \nonumber \nonumber \\
+ &&\exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \nonumber \nonumber \\
+ &&\left.\exp\left(iL_2 {\Delta t}/2\right)\right]^n \exp\left(iL_{\epsilon,2}(n{\Delta t}/2)\right) \nonumber \nonumber \\
+ &&\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right)\exp\left(iL_{\mathrm{NHC-baro}}(n{\Delta t}/2)\right),\end{aligned}
+ :label: eqnVVNpressure
where the box volume integration occurs every step, but the auxiliary
variable integrations happen every :math:`n` steps.
algorithm for the velocities becomes
.. math:: \mathbf{v}(t+{\frac{\Delta t}{2}})~=~\mathbf{f}_g * \lambda * \left[ \mathbf{v}(t-{\frac{\Delta t}{2}}) +\frac{\mathbf{F}(t)}{m}\Delta t + \mathbf{a}_h \Delta t \right],
+ :label: eqntotalupdate
where :math:`g` and :math:`h` are group indices which differ per atom.
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff