Микро- и макропараметры

Для описания МД систем применяют методы статистической механики, биологические системы обычно моделируют в NPT-ансамбле, в нём число частиц, давление и температура принимаются неизменными.

Температура

Мгновенную полную энергию системы \(E\) можно определить простым суммированием по кинетическим энергиям и по полным потенциалам взаимодействия двух частиц: \[ E(t)=\sum_{i}\dfrac{m_{i}{\bf v}_{i}^{2}(t)}{2}+\sum_{i<j}U_{i,j}, \] каждое парное взаимодействие входит в выражение лишь единожды.

Температура непосредственно связана со средней (по частицам) кинетической энергией, приходящейся на одну степень свободы: \[ T(t)=\dfrac{2}{3kN}\sum_{i}\dfrac{m_{i}{\bf {\bf v}}_{i}^{2}(t)}{2}, \tag{3.1}\] видно, что температура неявным образом зависит от времени, что нарушает предпосылки NPT-ансамбля, это противоречие можно устранить введя ограничение на среднюю кинетическую скорость.

Согласно Ур. 3.1, разность мгновенной \(T(t)\) и фиксированной в ансамбле \(\Theta\) температур

\[ \Delta T=T(t)-\Theta=\dfrac{2}{3kN}\left\{ \sum_{i}\lambda^{2}{\bf v}_{i}^{2}(t)-\sum_{i}{\bf v}_{i}^{2}(t)\right\}, \] коэффициент \(\lambda\) — некий параметр, показывающий среднее отклонение мгновенных скоростей от ожидаемых в условиях NPT-ансамбля, выделив явно: \[ \lambda=\sqrt{\dfrac{\Theta}{T(t)}}, \] этот параметр масштабирования скоростей включается на каждом шаге интегрирования: \[ {\bf v}_{i}(t+\delta t)\gets\lambda{\bf v}_{i}(t)+{\bf a}_{i}(t)\delta t\text{.} \]

Существуют и другие подходы поддержания постоянной температуры в системе, они воспроизводят теплообмен с внешним источником тепловой энергии — термостатом, а сами называются методами термостатирования.

Термостат Берендсена
[1] так же вводит масштабирующий коэффициент скорости \(\lambda_\t{B}\) для каждого шага интегрирования: \[\lambda_\t{B}=\sqrt{1+\dfrac{\delta t}{Q}\left(\dfrac{\Theta}{T(t)}-1\right)},\] интенсивность теплообмена между системой и термостатом учитывается выбором \(Q\).
Термостат Нозе-Гувера
[2,3] модифицирует уравнение движения Ур. 1.1, масштабируя координаты и скорости системы через параметр \[\lambda_\t{NG}=\dfrac{1}{Q}\left\{ \sum_{i}\dfrac{m_{i}{\bf v}_{i}^{2}(t)}{2}-3kN\Theta\right\} ,\] где \(Q\) учитывает интенсивность теплового обмена; задача МД сводится к решению систем уравнений движения вида \[\dv[2]{\bf r_{i}(t)}{t}=\dfrac{{\bf F_{i}}}{m_{i}}-\lambda_\t{NG}{\bf v_{i}},\] \(\lambda_\t{NG}\) по сути является коэффициентом вязкого трения.
Термостат Ланжевена
[46] основан на уравнениях движения броуновской динамики — вводятся две дополнительные силы: случайная составляющая \(\bf a\), вызывающая нагрев частиц, и диссипативная сила \(\bf b\), рассеивающая лишнюю энергию (по сути сила трения). Аналогично Ур. 1.1, уравнение движения: \[m_{i}\dv[2]{\bf r_{i}(t)}{t}=\sum_{i\ne j}{\bf F_{i,j}(t)+{\bf a_{i}+{\bf b_{i}\quad i,j=1..N,}}}\] случайная составляющая распределена нормально: \[\bf a_i \sim \cl N (\bf 0, \bm{\sigma}_i) \quad \t{и} \quad \sigma_i=\dfrac{m_i k \Theta}{\delta t \cdot Q},\] где вектора \(\mathbf{0},\bm{\sigma}_{i}\in\mathbb{R}^{3}\), \(\sigma_{i}\) — компонента \(\bm{\sigma}\), \(Q\) — параметр теплового обмена, \(\Theta\) — поддерживаемая температура. Для трения: \[\mathbf{b}_{i}=-Q\dfrac{m_{i}}{2}\dv{\bf r_{i}}{t}.\]

Давление

Для описания макропараметров реальных систем одни параметры раскладывают в степенной ряд по другим, например, по обратным степеням молярного объёма: \[ \dfrac{P}{RT}=\sum_{n=1}^{\infty}\dfrac{B_{i}(T)}{V_{m}^{n-1}}, \] где \(B_{i}(T)\) — функция, зависящая от температуры; такие уравнения называют вириальными. Вириальное уравнение не имеет чёткого физического смысла. Аналогичный подход используется и в МД, здесь мгновенное давление определяется вириальным соотношением

\[P(t)=\dfrac{1}{3V}\sum_{i}\left\{ \dfrac{{\bf p}_{i}^{2}(t)}{m_{i}}+\sum_{j\ne i}{\bf F}_{i}(t)\cdot{\bf r_{i,j}}\right\} ,\] что, как и в случае температуры, противоречит предпосылкам NPT-ансамбля. Эту проблему решают алгоритмы баростатирования.

Баростат Берендсена

[1] использует принцип локального возмущения, где изменение давления во времени определяется как \[\dv{P(t)}{t}=\dfrac{\Phi-P(t)}{Q},\] параметр \(Q\) учитывает интенсивность баростата, \(\Phi\) — поддерживаемое им давление путём изменения размера системы: объём масштабируется с помощью коэффициента \[\xi_{B}=\dfrac{\beta_{T}}{3Q}\left\{ P(t)-\Phi\right\} ,\] изменение объёма записывается как \[\Delta V=V(t)-V, \tag{3.2}\] уравнения движения принимают вид: \[\dv{\bf r_{i}(t)}{t}={\bf v}_{i}-\xi_{B}{\bf r_{i}}.\]

Баростат Нозе-Гувера

[2,3] аналогично вводит масштабирующий множитель: \[\dv{\bf \xi_{\text{NG}}}{t}=\dfrac{V(t)}{Q}\left\{ P(t)-\Phi\right\} ,\] изменение объёма как в Ур. 3.2, уравнения движения: \[\dv[2]{\bf r_{i}(t)}{t}=\dfrac{{\bf F_{i}}}{m_{i}}-\xi_{\text{NG}}{\bf v_{i}}.\]

Баростат Ланжевена

[7] основан на уравнениях Ланжевена и броуновской динамике. Принципы, лежащие в основе приведены в Глава 3.1, ввиду его сложных построений, ограничимся лишь указанием.

Размер

Если не наложены ограничения на размеры системы, она сама заполнит предоставленный ей объём, так можно описывать молекулы в вакууме или в газе, однако в ходе такого моделирования частицы разлетятся на огромные расстояния, и нам рано или поздно не хватит размера машинных чисел даже для простого задания координат; поэтому разумным приближением будет введение некой ячейки моделирования, подойдя к границе которой частицы испытают либо упругое соударение, либо их координаты будут транслированы на противоположную сторону этой ячейки.

Периодические граничных условия [8] могут задаваться прямоугольным параллелепипедом, гексагональной призмой, усечённым октаэдром, ромбическим додекаэдром; непериодические граничные условия можно задавать сферой.

В конденсированных фазах часто присутствует ближний и дальний порядок, это означает, что для описания всей системы достаточно описания небольшой её части. Здесь применение периодичности особенно полезно. Кроме того, за счёт изменения размеров ячейки моделирования возможно поддержание постоянного давления (Глава 3.2).

Источники

1.
Berendsen H.J.C. и др. Molecular dynamics with coupling to an external bath // The Journal of Chemical Physics. 1984. Т. 81, № 8. С. 3684–3690.
2.
Nosé S. A unified formulation of the constant temperature molecular dynamics methods // The Journal of Chemical Physics. 1984. Т. 81, № 1. С. 511–519.
3.
Hoover W.G. Canonical dynamics: Equilibrium phase-space distributions // Physical Review A. 1985. Т. 31, № 3. С. 1695–1697.
4.
Подрыга В.О., Поляков С.В. Молекулярно-динамическое моделирование процесса установления термодинамического равновесия нагретого никеля // Препринты Института прикладной математики им. М. В Келдыша РАН. 2014. С. 20–41.
5.
Davidchack R.L., Handel R., Tretyakov M.V. Langevin thermostat for rigid body dynamics // Journal of Chemical Physics. 2009. Т. 130, № 23.
6.
Brünger A.T. X-plor version 3.1: A system for x-ray crystallography and NMR // Нью-Хейвен (Коннектикут): Yale Univ Press. 1992. С. 39.
7.
Feller S.E. и др. Constant pressure molecular dynamics simulation: The Langevin piston method // The Journal of Chemical Physics. American Institute of Physics, 1995. Т. 103, № 11. С. 4613–4621.
8.
Bekker H., Van Den Berg J.P., Wassenaar T.A. A method to obtain a near-minimal-volume molecular simulation of a macromolecule, using periodic boundary conditions and rotational constraints // Journal of Computational Chemistry. 2004. Т. 25, № 8. С. 1037–1046.