Молекулярные взаимодействия в классической молекулярной динамике

1 Уравнение движения

В классической молекулярной динамике (МД) молекулы¹ представляются в виде идеализированных систем, состоящих из материальных точек — абстрактных моделей атомов и ионов. Эти частицы взаимодействуют друг с другом явно, т/е выражения соответствующих сил задаются явно. Полный набор аналитических соотношений вместе со значениями всех необходимых численных параметров, необходимых для расчёта сил, составляет силовое поле (СП), т/е силовое поле — это набор уравнений и констант, описывающих взаимодействия между молекулами.

¹ под молекулами подразумеваются втч атомы и ионы

Задачей МД является нахождение траекторий всех \(N\) частиц, входящих в систему, через решение 2го уравнения динамики Ньютона: \[ m_i \frac{d\mathbf{r}_i(t)}{dt} = \sum_{i \neq j}\mathbf{F}_{i,j}(t) \quad i,j=1..N, \tag{1}\] где \(m_i\) и \({\mathbf r_{i}}\) — масса и радиус-вектор \(i\)ой частицы (атома, иона идр), \({\mathbf F}_{i,j}\) — сила действующая на \(i\)ю частицу со стороны \(j\)ой частицы, величину \({\mathbf F}_{i,j}\) определяет потенциал: \[ \mathbf{F}_{i,j} = -\nabla U_{i,j}, \tag{2}\] где \(U_{i,j}\) — потенциал взаимодействия частиц, \(\nabla\) — оператор набла, \(\nabla U_{i,j}\) — градиент \(U_{i,j}\).

Полная потенциальная энергия взаимодействия частиц \(U\) складывается из нескольких составляющих: \[ U = \underbrace{U_b + U_v + U_t + U_f}_{\text{ковалентные}} + \underbrace{U_q + U_d}_{\text{нековалентные}}. \tag{3}\]

Далее приводятся конкретные выражения термов из 3, используемые в большинстве подходах МД. В разных СП эти термы могут вычисляться по-разному, однако, различия скорее касаются деталей.

flowchart LR
    I[взаимодействия]
        I --> C[ковалентные]
            C --> P[параболические]
                P --> L(длины связей)
                P --> A(валентные углы)
                P -..-> D(двугранные углы)
            C --> S[ряды]
                S --> D
        I --> N[нековалентные]
            N --> E(электростатические)
            N --> W[вандерваальсовы]
                W --> O(ориентационные)
                W --> Di(дисперсионные)
                W --> In(индукционные)

click L "#sec-bond-lengths"
click A "#sec-valence-angles"
click D "#sec-dihedrals"
click E "#sec-electrostatics"
click W "#sec-vanderwaals"

2 Ковалентные взаимодействия

Параболические

Потенциалы длин связей 4 и валентных углов 5 приближаются законом Гука (параболой), хотя их природе соответствует потенциал Ленарда—Джонса 10. Упрощённые выражения используются из-за простоты их численного вычисления и другой важной особенности: при удалении длины связи \(l\) от равновесного значения \(\bar l\) на бесконечность, стремительно растёт и соответствующий параболический потенциал \(U_b \propto (l - \bar l)^2\). В этом приближении смоделировать разрыв связи невозможно в принципе, и это хорошо — в процессе численного эксперимента целостность молекулы не нарушится. Такой системе не страшны даже высокие температуры.

Длины связей

Потенциал — аддитивная величина. Суммарная энергия деформации длин ковалентных связей: \[U_b = \frac{1}{2} \sum_i K_{b,i} (l_i - \bar{l}_i)^2, \tag{4}\] где \(i\) — номер связи, \(K_{b,i}\) — эффективная жёсткость гуковской «пружины», \(l_i\) и \(\bar l_i\) — длина и равновесная длина связи;

Валентные углы

\[U_v = \frac{1}{2}\sum_i K_{v,i}(\alpha_i - \bar{\alpha}_i)^2, \tag{5}\] где \(K_{v,i}\) — эффективная упругость валентного угла \(\alpha_i\), \(\bar \alpha_i\) — его равновесное значение;

Ряды

Торсионные (двугранные) углы

Потенциалы двугранных углов 6 и плоских групп 7 представляются рядами Фурье, как правило, для корректного описания достаточно 4х членов ряда (включая нулевой)

Периодические потенциалы вращения вокруг простых связей описываемых разложением в ряд Фурье \[ U_t = \frac{1}{2} \sum_i \sum_l K_{\phi,l} \left(1 + g_{i,l} \cos(n_{i,l}\phi_i) \right), \tag{6}\] суммирование ведётся по номерам торсионных углов \(i\) и гармоникам \(l\), константа \(K_{\phi,l}\) — аналог жёсткости, \(g_{i,l} \in [-1\,..\,1]\) — вклад гармоники, \(n_{i,l}\) — её кратность.

Плоские группы

Для описания плоских фрагментов, например, пептидных связей, достаточно задать необходимое количество двухгранных углов между атомами, от которых требуется находиться в одной плоскости, и приписать соответствующие потенциалы с высоким барьером вращения, препятствующим выходу атомов из плоскости: \[U_{f} = \cdots, \tag{7}\] задаётся аналогично 6, с поправкой на собственную константу \(K_{f,l}\);

3 Нековалентные взаимодействия

кулоновские 8, вандерваальсовы 10 и водородные 11 взаимодействия суммируются по принципу «каждый с каждым»

Кулоновские (электростатические) взаимодействия

\[ U_q = \frac{1}{4\pi\epsilon\epsilon_0}\sum_{i < j}\frac{q_i q_j}{r_{i,j}}, \tag{8}\]

здесь и ниже \(r_{i,j}=|{\mathbf r_{i}-{\mathbf r_{j}|}}\), \(q_{i}\) и \(q_{j}\) — заряды частиц, \(\epsilon\) и \(\epsilon_{0}\) — диэлектрические проницаемости среды и вакуума, \(i < j\) подчёркивает, что взаимодействие между двумя частицами считается лишь единожды;

Кулоновский потенциал 8 медленно убывает на расстоянии (\(U_{q}\propto r^{-1}\)), поэтому нередко применяются 2 подхода для уменьшения учёта этого взаимодействия:

• вводится радиус отсечения — расстояние за границами которого от попарного суммирования электростатических взаимодействий переходят к суммированию по Эвальду Yeh and Berkowitz (1999) (обычно применяется для периодических систем в конденсированном состоянии), а порой вовсе пренебрегают этими дальнодействующими взаимодействиями;
• диэлектрическая проницаемость вакуума \(\epsilon_{0}\) в 8 заменяется на расстояние \(r\), в результате потенциал спадает быстрее (\(U_{q}\propto r^{-2}\)).
Yeh, In Chul, and Max L. Berkowitz. 1999. “Ewald Summation for Systems with Slab Geometry.” Journal of Chemical Physics 111 (7): 3155–62. https://doi.org/10.1063/1.479595.

Радиус отсечения может вводиться и для иных типов взаимодействий, суммирующихся по принципу «каждый с каждым», всё потому, что асимптотическая сложность подобного суммирования \(O(N^{2})\)

Вандерваальсовы взаимодействия

Простейший диполь состоит из двух материальных точек с одинаковыми по величине разноимёными зарядами (\(+q\) и \(-q\)), находящихся на расстоянии \(l\) друг от друга. Для такой системы электрическим² дипольным моментом называют вектор³ \[ \mathbf p = q \mathbf l, \] где \(\mathbf l\) — вектор с началом в \(-q\) и концом в \(+q\)⁴.

² ёщё бывает магнитный

³ электрический дипольный момент молекул традиционно измеряют в единицах системы СГС — в дебаях (\(\mathrm Д\), \(\mathrm D\)), эквивалентной единицей в СИ является метр \(\times\) кулон

⁴ дипольный момент направлен от \(\ominus\) к \(\oplus\), в отличие от линий напряжённости \(\mathbf E\)

Дипольным моментом обладает любая система зарядов. В общем случае: \[ \mathbf p = \sum_i q_i \mathbf r_i, \] где \(q_i\) — величина \(i\)го заряда, \(\mathbf r_i\) — его положение.

Внешнее электрическое поле может приводить к перераспределению заряда в системе, что повлечёт изменение дипольных моментов. Величину пропорциональности между индуцированным (наведённым) дипольным моментом \(\mathbf p\) и напряжённостью внешнего поля \(\mathbf E\) называют поляризуемостью \(\alpha\): \[ \mathbf p = \alpha \mathbf E. \]

Силы межмолекулярного взаимодействия, возникающие в результате взаимодействия диполей (истинных и наведённых) называют силами Ван-дер-Ваальса⁵, их разделяются на 3 типа:

⁵ Взяв уравнения состояния идеального газа Ван-дер-Ваальс предложил уравнение для реального газа, добавив 2 параметра (\(a\) и \(b\)) для учёта притяжения и отталкивания между молекулами: \[ \biggl( P + \underbrace{\frac{a}{V_m^2}}_{P_i} \biggr) \left(V_m - b\right) = RT. \]

• \(b\) — объём, занимаемый молекулами газа и поэтому недоступный из-за взаимного отталкивания.
• \(a\) — параметр притяжения молекул, который проявляется в возникновении избыточного давлении \(P_i\): силы межмолекулярного притяжения стремятся собрать газ в меньший объём, они действуют из толщи газа и отсутствуют на стенках. Слагаемое \(P\) обусловленно столкновениями.

1. ориентационные (диполь-дипольные) — взаимодействие между двумя истинными диполями;
2. индукционные (поляризационные, дебаевские) — между истинным и наведённым диполем;
3. дисперсионные (лондоновские) — между двумя наведёнными диполями.

В классической МД ориентационную составляющую не требуется описывать специально, поскольку для расчёта сил взаимодействия двух диполей достаточно закона Кулона 8, выражения для потенциала пружины 4 между положительным и отрицательным концами диполя и 2го закона Ньютона 1.

В первом приближении поляризуемостью атомов можно принебречь и не учитывать поляризационные взаимодействия. Биополимеры в основном состоят из слабополяризуемых атомов 1го⁶ и 2го периодов: \(\text{H}\), \(\text{C}\), \(\text{N}\), \(\text{O}\); заметной⁷ поляризуемостью обладают атомы халькогенов (\(\text{S}\), \(\text{Se}\)), тяжёлых галогенов (\(\text{Br}\), \(\text{I}\)), \(d\)- и \(f\)-элементов, иногда соответствующие им катионы (\(\text{Cs}^+\), \(\text{Rb}^+\), \(\text{Ba}^{2+}\)). Анионы обладают большей поляризуемостью (\(\text{O}^-\), \(\text{S}^-\), \(\text{Br}^-\)).

⁶ атом \(\text{H}\) может быть как неполяризован (в связи \(\text{CH}\)), так и поляризован (\(\text{OH}\), \(\text{NH}\) итд)

⁷ если поляризуемостью нельзя принебречь, следует использовать специальные силовые поля

Индукционное взаимодействие обусловлено мгновенными дипольными моментами, самопроизвольно возникающими у атомов и молекул в следствие флуктуаций электронной плотности и перераспределения заряда. Случайно возникший диполь индуцирует перераспределение электронной плотности в соседних атомах, молекулах и ионах, приводя к появлению у них собственного дипольного момента. Этот тип взаимодействия универсален, т/е проявляется во всех случаях, поскольку возникает самопроизвольно между любыми близкими частицами. В случае ионов, однако, его сила будет малозаметна на фоне более сильного притяжения или отталкивания истинных зарядов.

Индукционные взаимодействия хотя и являются слабыми, играют большую роль: именно они отвечают за гидрофобный эффект, возникающий при контакте 2х фаз (полярной и неполярной), что невозможно не учитывать при моделировании белковых глобул, липидных слоёв и, например, там, где индукционное взаимодействие будет основным или даже единственным — в сжиженных углеводородах или в благородных газах.

Таким образом, вандерваальсовы силы часто сводят к индукционным и приближают их потенциалом Ленарда — Джонса: \[ U_{LJ}(r) = 4\epsilon\left\{\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right\}, \tag{9}\] где \(\epsilon\) — глубина потенциальной ямы, \(\sigma\) — расстояние, при котором \(U_{LJ}(\sigma)=0\); степенная функция \(r^{-6}\) аппроксимирует притяжение частиц, \(r^{-12}\) — отталкивание, отсюда это выражение также называют «потенциалом 6-12», он описывает энергию взаимодействия частиц как функцию расстояния между ними.

Существуют и другие параметризации 9: \[ U_d = \sum_{i < j} \left\{\frac{A_{i,j}}{r_{i,j}^{12}} - \frac{B_{i,j}}{r_{i,j}^6}\right\}, \tag{10}\] также иногда используется потенциал Букиннгема: \[ U_d = \sum_{i < j} \left \{ A_{i,j} \exp (- B_{i,j} r_{i,j}) - \dfrac{C_{i,j}}{r_{i,j}^6} \right \}, \] где \(A_{i,j}\), \(B_{i,j}\) и \(C_{i,j}\) — константы, параметризующие потенциальную яму для пары взаимодействующих частиц \(i\) и \(j\).

Водородная связь

Ввиду своей специфичности, вызванной малыми размерами и высокой плотностью заряда на атоме водорода, иногда полярному водороду приписывают потенциал \[ U_{h}= \sum_{i < j} \left \{ \dfrac{A_{i,j}'}{r_{i,j}^{12}}-\dfrac{B_{i,j}'}{r_{i,j}^{10}} \right \}, \tag{11}\] \(A_{i,j}'\) и \(B_{i,j}'\) — аналогично 10 с поправкой на степень при \(B_{i,j}'\). Член \(r^{-10}\) берётся из феноменологических соображений.