Оценка термодинамических микро- и макропараметров в классической молекулярной динамике

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

1 Температура

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

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

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

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

Термостат Берендсена (Berendsen et al. 1984)

так же вводит масштабирующий коэффициент скорости \(\lambda_\text{B}\) для каждого шага интегрирования: \[\lambda_\text{B} = \sqrt{1 + \frac{\delta t}{Q} \left( \frac{\Theta}{T(t)} - 1 \right)},\] интенсивность теплообмена между системой и термостатом учитывается выбором \(Q\).

Термостат Нозе-Гувера (Nosé 1984; Hoover 1985)

модифицирует уравнение движения, масштабируя координаты и скорости системы через параметр \[\lambda_\text{NG} = \frac{1}{Q} \left\{ \sum_{i} \frac{m_{i} \mathbf{v}_{i}^{2}(t)}{2} - 3kN\Theta \right\},\] где \(Q\) учитывает интенсивность теплового обмена; задача МД сводится к решению систем уравнений движения вида \[\frac{d^2 \mathbf{r}_{i}(t)}{dt^2} = \frac{\mathbf{F}_{i}}{m_{i}} - \lambda_\text{NG} \mathbf{v}_{i},\] \(\lambda_\text{NG}\) по сути является коэффициентом вязкого трения.

Термостат Ланжевена (Подрыга and Поляков 2014; Davidchack, Handel, and Tretyakov 2009; Brünger 1992)

Подрыга, Виктория Олеговна, and Сергей Владимирович Поляков. 2014. “Молекулярно-Динамическое Моделирование Процесса Установления Термодинамического Равновесия Нагретого Никеля.” Препринты Института Прикладной Математики Им. М. В Келдыша РАН, 20–41.

Davidchack, Ruslan L., Richard Handel, and M. V. Tretyakov. 2009. “Langevin Thermostat for Rigid Body Dynamics.” Journal of Chemical Physics 130 (23). https://doi.org/10.1063/1.3149788.

Brünger, Axel T. 1992. “X-Plor Version 3.1: A System for x-Ray Crystallography and NMR.” Нью-Хейвен (Коннектикут): Yale Univ Press, 39.

основан на уравнениях движения броуновской динамики — вводятся две дополнительные силы: случайная составляющая \(\mathbf{a}\), вызывающая нагрев частиц, и диссипативная сила \(\mathbf{b}\), рассеивающая лишнюю энергию (по сути сила трения). Аналогично уравнению Ньютона, уравнение движения: \[m_{i} \frac{d^2 \mathbf{r}_{i}(t)}{dt^2} = \sum_{i \ne j} \mathbf{F}_{i,j}(t) + \mathbf{a}_{i} + \mathbf{b}_{i} \quad i,j=1..N,\] случайная составляющая распределена нормально: \[\mathbf{a}_{i} \sim \mathcal{N} (\mathbf{0}, \boldsymbol{\sigma}_{i}) \quad \text{и} \quad \sigma_{i} = \frac{m_{i} k \Theta}{\delta t \cdot Q},\] где вектора \(\mathbf{0}, \boldsymbol{\sigma}_{i} \in \mathbb{R}^{3}\), \(\sigma_{i}\) — компонента \(\boldsymbol{\sigma}\), \(Q\) — параметр теплового обмена, \(\Theta\) — поддерживаемая температура. Для трения: \[\mathbf{b}_{i} = -Q \frac{m_{i}}{2} \frac{d \mathbf{r}_{i}}{dt}.\]

2 Давление

Для описания макропараметров реальных систем одни параметры раскладывают в степенной ряд по другим, например, по обратным степеням молярного объёма: \[ \frac{P}{RT} = \sum_{n=1}^{\infty} \frac{B_{i}(T)}{V_{m}^{n-1}}, \] где \(B_{i}(T)\) — функция, зависящая от температуры; такие уравнения называют вириальными. Вириальное уравнение не имеет чёткого физического смысла. Аналогичный подход используется и в МД, здесь мгновенное давление определяется вириальным соотношением \[ P(t) = \frac{1}{3V} \sum_{i} \left\{ \frac{\mathbf{p}_{i}^{2}(t)}{m_{i}} + \sum_{j \ne i} \mathbf{F}_{i}(t) \cdot \mathbf{r}_{i,j} \right\}, \] что, как и в случае температуры, противоречит предпосылкам NPT-ансамбля. Эту проблему решают алгоритмы баростатирования.

Баростат Берендсена (Berendsen et al. 1984)

Berendsen, H. J.C., J. P.M. Postma, W. F. Van Gunsteren, A. Dinola, and J. R. Haak. 1984. “Molecular Dynamics with Coupling to an External Bath.” The Journal of Chemical Physics 81 (8): 3684–90. https://doi.org/10.1063/1.448118.

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

Баростат Нозе-Гувера (Nosé 1984; Hoover 1985)

Nosé, Shuichi. 1984. “A Unified Formulation of the Constant Temperature Molecular Dynamics Methods.” The Journal of Chemical Physics 81 (1): 511–19. https://doi.org/10.1063/1.447334.

Hoover, William G. 1985. “Canonical Dynamics: Equilibrium Phase-Space Distributions.” Physical Review A 31 (3): 1695–97. https://doi.org/10.1103/PhysRevA.31.1695.

аналогично вводит масштабирующий множитель: \[\frac{d \xi_{\text{NG}}}{dt} = \frac{V(t)}{Q} \left\{ P(t) - \Phi \right\},\] изменение объёма как в 2, уравнения движения: \[\frac{d^2 \mathbf{r}_{i}(t)}{dt^2} = \frac{\mathbf{F}_{i}}{m_{i}} - \xi_{\text{NG}} \mathbf{v}_{i}.\]

Баростат Ланжевена (Feller et al. 1995)

Feller, Scott E., Yuhong Zhang, Richard W. Pastor, and Bernard R. Brooks. 1995. “Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method.” The Journal of Chemical Physics 103 (11): 4613–21. https://doi.org/10.1063/1.470648.

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

3 Размер

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

Периодические граничные условия (Bekker, Van Den Berg, and Wassenaar 2004) могут задаваться прямоугольным параллелепипедом, гексагональной призмой, усечённым октаэдром, ромбическим додекаэдром; непериодические граничные условия можно задавать сферой.

Bekker, Henk, Jur P. Van Den Berg, and Tsjerk A. Wassenaar. 2004. “A Method to Obtain a Near-Minimal-Volume Molecular Simulation of a Macromolecule, Using Periodic Boundary Conditions and Rotational Constraints.” Journal of Computational Chemistry 25 (8): 1037–46. https://doi.org/10.1002/jcc.20050.

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