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

molecular dynamics
russian
Published

January 1, 2021

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

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

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

  • вводится радиус отсечения — расстояние за границами которого от попарного суммирования электростатических взаимодействий переходят к суммированию по Эвальду [1] (обычно применяется для периодических систем в конденсированном состоянии), а порой вовсе пренебрегают этими дальнодействующими взаимодействиями;
  • диэлектрическая проницаемость вакуума \epsilon_{0} в 8 заменяется на расстояние r, в результате потенциал спадает быстрее (U_{q}\propto r^{-2}).
    • 1.
    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 друг от друга. Для такой системы электрическим2 дипольным моментом называют вектор3 \mathbf p = q \mathbf l, где \mathbf l — вектор с началом в -q и концом в +q4.

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

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

4 дипольный момент направлен от \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.

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

5 Взяв уравнения состояния идеального газа Ван-дер-Ваальс предложил уравнение для реального газа, добавив 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го6 и 2го периодов: \text{H}, \text{C}, \text{N}, \text{O}; заметной7 поляризуемостью обладают атомы халькогенов (\text{S}, \text{Se}), тяжёлых галогенов (\text{Br}, \text{I}), d- и f-элементов, иногда соответствующие им катионы (\text{Cs}^+, \text{Rb}^+, \text{Ba}^{2+}). Анионы обладают большей поляризуемостью (\text{O}^-, \text{S}^-, \text{Br}^-).

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

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

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

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