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"
Молекулярные взаимодействия
Уравнение движения
В классической молекулярной динамике (МД) молекулы1 представляются в виде идеализированных систем, состоящих из материальных точек — абстрактных моделей атомов и ионов. Эти частицы взаимодействуют друг с другом явно, т/е выражения соответствующих сил задаются явно. Полный набор аналитических соотношений вместе со значениями всех необходимых численных параметров, необходимых для расчёта сил, составляет силовое поле (СП), т/е силовое поле — это набор уравнений и констант, описывающих взаимодействия между молекулами.
1 под молекулами подразумеваются втч атомы и ионы
Задачей МД является нахождение траекторий всех \(N\) частиц, входящих в систему, через решение 2го уравнения динамики Ньютона: \[ m_{i}\dv[2]{\bf r_{i}(t)}{t}=\sum_{i \ne j}{\bf F_{i,j}(t) \qquad i,j=1..N,} \tag{1.1}\] где \(m_i\) и \({\bf r_{i}}\) — масса и радиус-вектор \(i\)ой частицы (атома, иона идр), \({\bf F}_{i,j}\) — сила действующая на \(i\)ю частицу со стороны \(j\)ой частицы, величину \({\bf F}_{i,j}\) определяет потенциал: \[ \mathbf{F}_{i,j} = -\nabla U_{i,j}, \tag{1.2}\] где \(U_{i,j}\) — потенциал взаимодействия частиц, \(\nabla\) — оператор набла, \(\nabla U_{i,j}\) — градиент \(U_{i,j}\).
Полная потенциальная энергия взаимодействия частиц \(U\) складывается из нескольких составляющих: \[ U = \ub{U_{b}+U_{v}+U_{t}+U_{f}}_\t{ковалентные}+\ub{U_{q}+U_{d}}_\t{нековалентные}. \tag{1.3}\]
Далее приводятся конкретные выражения термов из Ур. 1.3, используемые в большинстве подходах МД. В разных СП эти термы могут вычисляться по-разному, однако, различия скорее касаются деталей.
Ковалентные взаимодействия
Параболические
Потенциалы длин связей Ур. 1.4 и валентных углов Ур. 1.5 приближаются законом Гука (параболой), хотя их природе соответствует потенциал Ленарда—Джонса Ур. 1.10. Упрощённые выражения используются из-за простоты их численного вычисления и другой важной особенности: при удалении длины связи \(l\) от равновесного значения \(\bar l\) на бесконечность, стремительно растёт и соответствующий параболический потенциал \(U_b \propto (l - \bar l)^2\). В этом приближении смоделировать разрыв связи невозможно в принципе, и это хорошо — в процессе численного эксперимента целостность молекулы не нарушится. Такой системе не страшны даже высокие температуры.
Длины связей
Потенциал — аддитивная величина. Суммарная энергия деформации длин ковалентных связей: \[U_b = \dfrac 1 2 \sum_i K_{b,i} (l_{i} - \bar l_i)^2, \tag{1.4}\] где \(i\) — номер связи, \(K_{b,i}\) — эффективная жёсткость гуковской «пружины», \(l_i\) и \(\bar l_i\) — длина и равновесная длина связи;
Валентные углы
\[U_{v}=\dfrac{1}{2}\sum_{i}K_{v,i}(\alpha_{i}-\bar \alpha_i)^{2}, \tag{1.5}\] где \(K_{v,i}\) — эффективная упругость валентного угла \(\alpha_i\), \(\bar \alpha_i\) — его равновесное значение;
Ряды
Торсионные (двугранные) углы
Периодические потенциалы вращения вокруг простых связей описываемых разложением в ряд Фурье \[ U_t = \dfrac 1 2 \sum_i \sum_l K_{\phi,l} \left(1 + g_{i,l} \cos \{ n_{i,l}\phi_{i} \} \right), \tag{1.6}\] суммирование ведётся по номерам торсионных углов \(i\) и гармоникам \(l\), константа \(K_{\phi,l}\) — аналог жёсткости, \(g_{i,l} \in [-1\,..\,1]\) — вклад гармоники, \(n_{i,l}\) — её кратность.
Плоские группы
Для описания плоских фрагментов, например, пептидных связей, достаточно задать необходимое количество двухгранных углов между атомами, от которых требуется находиться в одной плоскости, и приписать соответствующие потенциалы с высоким барьером вращения, препятствующим выходу атомов из плоскости: \[ U_{f} = \cdots, \tag{1.7}\] задаётся аналогично Ур. 1.6, с поправкой на собственную константу \(K_{f,l}\);
Нековалентные взаимодействия
Кулоновские (электростатические) взаимодействия
\[ U_{q}=\dfrac{1}{4\pi\epsilon\epsilon_{0}}\sum_{i < j}\dfrac{q_{i}q_{j}}{r_{i,j}}, \tag{1.8}\] здесь и ниже \(r_{i,j}=|{\bf r_{i}-{\bf r_{j}|}}\), \(q_{i}\) и \(q_{j}\) — заряды частиц, \(\epsilon\) и \(\epsilon_{0}\) — диэлектрические проницаемости среды и вакуума, \(i < j\) подчёркивает, что взаимодействие между двумя частицами считается лишь единожды;
Кулоновский потенциал Ур. 1.8 медленно убывает на расстоянии (\(U_{q}\propto r^{-1}\)), поэтому нередко применяются 2 подхода для уменьшения учёта этого взаимодействия:
- вводится радиус отсечения — расстояние за границами которого от попарного суммирования электростатических взаимодействий переходят к суммированию по Эвальду [1] (обычно применяется для периодических систем в конденсированном состоянии), а порой вовсе пренебрегают этими дальнодействующими взаимодействиями;
- диэлектрическая проницаемость вакуума \(\epsilon_{0}\) в Ур. 1.8 заменяется на расстояние \(r\), в результате потенциал спадает быстрее (\(U_{q}\propto r^{-2}\)).
Радиус отсечения может вводиться и для иных типов взаимодействий, суммирующихся по принципу «каждый с каждым», всё потому, что асимптотическая сложность подобного суммирования \(O(N^{2})\)
Вандерваальсовы взаимодействия
Простейший диполь состоит из двух материальных точек с одинаковыми по величине разноимёными зарядами (\(+q\) и \(-q\)), находящихся на расстоянии \(l\) друг от друга. Для такой системы электрическим2 дипольным моментом называют вектор3 \[ \bf p = q \bf l, \] где \(\bf l\) — вектор с началом в \(-q\) и концом в \(+q\)4.
2 ёщё бывает магнитный
3 электрический дипольный момент молекул традиционно измеряют в единицах системы СГС — в дебаях (\(\u Д\), \(\u D\)), эквивалентной единицей в СИ является метр \(\times\) кулон
4 дипольный момент направлен от \(\ominus\) к \(\oplus\), в отличие от линий напряжённости \(\bf E\)
Код
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from sympy.abc import x, y
# electric potential of a point charge
def phi(x, y):
return 1 / sp.sqrt(x**2 + y**2)
# total potential
phi_tot = phi(x - 1/3, y) - phi(x + 1/3, y)
# electric field projections
E_x = -phi_tot.diff(x)
E_y = -phi_tot.diff(y)
# convert symbolic expressions to functions
E_x = sp.lambdify((x, y), E_x)
E_y = sp.lambdify((x, y), E_y)
phi_tot = sp.lambdify((x, y), phi_tot)
X, Y = np.meshgrid(np.linspace(-1, 1, 100), np.linspace(-1, 1, 100))
U = E_x(X, Y)
V = E_y(X, Y)
color = phi_tot(X, Y)
plt.imshow(
color,
extent=[-1, 1, -1, 1],
cmap='bwr',
interpolation='bicubic',
)
plt.streamplot(X, Y, U, V, density=1.5, arrowsize=0.5, linewidth=0.8)
plt.show()
Дипольным моментом обладает любая система зарядов. В общем случае: \[ \bf p = \sum_i q_i \bf r_i, \] где \(q_i\) — величина \(i\)го заряда, \(\bf r_i\) — его положение.
Внешнее электрическое поле может приводить к перераспределению заряда в системе, что повлечёт изменение дипольных моментов. Величину пропорциональности между индуцированным (наведённым) дипольным моментом \(\bf p\) и напряжённостью внешнего поля \(\bf E\) называют поляризуемостью \(\alpha\): \[ \bf p = \alpha \bf E. \]
Силы межмолекулярного взаимодействия, возникающие в результате взаимодействия диполей (истинных и наведённых) называют силами Ван-дер-Ваальса5, их разделяются на 3 типа:
5 Взяв уравнения состояния идеального газа Ван-дер-Ваальс предложил уравнение для реального газа, добавив 2 параметра (\(a\) и \(b\)) для учёта притяжения и отталкивания между молекулами: \[ \biggl( P + \ub{\frac{a}{V_m^2}}_{P_i} \biggr) \left(V_m - b\right) = RT. \]
- \(b\) — объём, занимаемый молекулами газа и поэтому недоступный из-за взаимного отталкивания.
- \(a\) — параметр притяжения молекул, который проявляется в возникновении избыточного давлении \(P_i\): силы межмолекулярного притяжения стремятся собрать газ в меньший объём, они действуют из толщи газа и отсутствуют на стенках. Слагаемое \(P\) обусловленно столкновениями.
- ориентационные (диполь-дипольные) — взаимодействие между двумя истинными диполями;
- индукционные (поляризационные, дебаевские) — между истинным и наведённым диполем;
- дисперсионные (лондоновские) — между двумя наведёнными диполями.
В классической МД ориентационную составляющую не требуется описывать специально, поскольку для расчёта сил взаимодействия двух диполей достаточно закона Кулона Ур. 1.8, выражения для потенциала пружины Ур. 1.4 между положительным и отрицательным концами диполя и 2го закона Ньютона Ур. 1.1.
В первом приближении поляризуемостью атомов можно принебречь и не учитывать поляризационные взаимодействия. Биополимеры в основном состоят из слабополяризуемых атомов 1го6 и 2го периодов: H, C, N, O; заметной7 поляризуемостью обладают атомы халькогенов (S, Se), тяжёлых галогенов (Br, I), \(d\)- и \(f\)-элементов, иногда соответствующие им катионы (Cs+, Rb+, Ba2+). Анионы обладают большей поляризуемостью (O−, S−, Br−).
6 атом H может быть как неполяризован (в связи CH), так и поляризован (OH, NH итд)
7 если поляризуемостью нельзя принебречь, следует использовать специальные силовые поля (Глава 4.3.1)
Индукционное взаимодействие обусловлено мгновенными дипольными моментами, самопроизвольно возникающими у атомов и молекул в следствие флуктуаций электронной плотности и перераспределения заряда. Случайно возникший диполь индуцирует перераспределение электронной плотности в соседних атомах, молекулах и ионах, приводя к появлению у них собственного дипольного момента. Этот тип взаимодействия универсален, т/е проявляется во всех случаях, поскольку возникает самопроизвольно между любыми близкими частицами. В случае ионов, однако, его сила будет малозаметна на фоне более сильного притяжения или отталкивания истинных зарядов.
Индукционные взаимодействия хотя и являются слабыми, играют большую роль: именно они отвечают за гидрофобный эффект, возникающий при контакте 2х фаз (полярной и неполярной), что невозможно не учитывать при моделировании белковых глобул, липидных слоёв и, например, там, где индукционное взаимодействие будет основным или даже единственным — в сжиженных углеводородах или в благородных газах.
Таким образом, вандерваальсовы силы часто сводят к индукционным и приближают их потенциалом Ленарда — Джонса: \[ U_{LJ}(r)=4\epsilon\left\{ \left(\dfrac{\sigma}{r}\right)^{12}-\left(\dfrac{\sigma}{r}\right)^{6}\right\}, \tag{1.9}\] где \(\epsilon\) — глубина потенциальной ямы, \(\sigma\) — расстояние, при котором \(U_{LJ}(\sigma)=0\); степенная функция \(r^{-6}\) аппроксимирует притяжение частиц, \(r^{-12}\) — отталкивание, отсюда это выражение также называют «потенциалом 6-12», он описывает энергию взаимодействия частиц как функцию расстояния между ними.
Существуют и другие параметризации Ур. 1.9: \[ U_{d}= \sum_{i < j} \left \{ \dfrac{A_{i,j}}{r_{i,j}^{12}}-\dfrac{B_{i,j}}{r_{i,j}^{6}} \right \}, \tag{1.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{1.11}\] \(A_{i,j}'\) и \(B_{i,j}'\) — аналогично Ур. 1.10 с поправкой на степень при \(B_{i,j}'\). Член \(r^{-10}\) берётся из феноменологических соображений.
Источники
1.
Yeh I.C., Berkowitz M.L. Ewald summation for systems with slab geometry // Journal of Chemical Physics. American Institute of Physics, 1999. Т. 111, № 7. С. 3155–3162.