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

Published

January 1, 2021

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, используемые в большинстве подходах МД. В разных СП эти термы могут вычисляться по-разному, однако, различия скорее касаются деталей.

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 подхода для уменьшения учёта этого взаимодействия:

вводится радиус отсечения — расстояние за границами которого от попарного суммирования электростатических взаимодействий переходят к суммированию по Эвальду [1] (обычно применяется для периодических систем в конденсированном состоянии), а порой вовсе пренебрегают этими дальнодействующими взаимодействиями;
диэлектрическая проницаемость вакуума $\epsilon_{0}$ в 8 заменяется на расстояние $r$ , в результате потенциал спадает быстрее ( $U_{q}\propto r^{-2}$ ).

Yeh I.C., Berkowitz M.L. Ewald summation for systems with slab geometry // Journal of Chemical Physics. American Institute of Physics, 1999. Vol. 111, № 7. P. 3155–3162.

Радиус отсечения может вводиться и для иных типов взаимодействий, суммирующихся по принципу «каждый с каждым», всё потому, что асимптотическая сложность подобного суммирования $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$ обусловленно столкновениями.

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

В классической МД ориентационную составляющую не требуется описывать специально, поскольку для расчёта сил взаимодействия двух диполей достаточно закона Кулона 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}$ берётся из феноменологических соображений.