Интегрирование уравнений движения в молекулярной динамике

1 Алгоритм Верле

Задача МД сводится к решению дифференциального уравнения Ньютона для каждой частицы, всего таких уравнений \(N\) в векторной форме или \(3N\) в скалярной (для трёхмерного пространства), решением системы будут \(3N\) проекций координат¹ \(r_{i,x}(t)\), \(r_{i,y}(t)\), \(r_{i,z}(t)\) и \(3N\) проекций скоростей \(v_{i,x}(t)\), \(v_{i,y}(t)\), \(v_{i,z}(t)\) для всех частиц системы в каждый момент времени.

¹ в данном случае декартовых

Verlet, Loup. 1967. “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules.” Physical Review 159 (1): 98–103. https://doi.org/10.1103/PhysRev.159.98.

Известно, что аналитический аппарат классической механики не может справиться с решением подобной задачи уже при \(N=3\) (задача трёх тел), поэтому здесь применимы только численные подходы: временна́я переменная \(t\) принимается дискретной, прирастающей на шаг интегрирования \(\delta t\), от решения дифференциальных уравнений переходят к их численному интегрированию, для этого обычно применяется алгоритм Верле Verlet (1967) и его модификации; в МД необходимо знать не только положение частиц, но и их скорости; последние можно определить, если использовать этот алгоритм в т.н. скоростной форме.

\begin{algorithm} \caption{Интегрирование методом Верле} \begin{algorithmic} \For{ частицы $i$ с 1ой по $N$ую} \State $\mathbf r_i(t + \delta t) \gets \mathbf r_i(t) + \mathbf v_i(t) \delta t + \mathbf a_i(t) \dfrac{\delta t^2}{2}$\Comment{$\mathbf r_i$, $\mathbf v_i$, $\mathbf a_i$ c предыдущего шага} \State $\mathbf v_i(t + \delta t / 2) \gets \mathbf v_i(t) + \mathbf a_i(t) \dfrac{\delta t}{2}$\Comment{скорость на половинном шаге} \State $ \mathbf a_i(t + \delta t) \gets \dfrac{\mathbf F_i (t + \delta t)}{m_i}$\Comment{$\displaystyle{\mathbf F_i = -\sum_{j \ne i} \frac{\partial U(r_{i,j})}{\partial r_{i,j}} \dfrac{\mathbf r_{i,j}}{r_{i,j}}}$} \State $\mathbf v_i(t + \delta t) \gets \mathbf v_i(t + \delta t/2) + \mathbf a_i(t + \delta t) \dfrac{\delta t}{2}$\Comment{скорость на целом шаге} \EndFor \end{algorithmic} \end{algorithm}

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

2 Фазовое пространство

Результатом интегрирования системы уравнений Ньютона будет набор скоростей и координат всех частиц системы, зависящих от времени; их можно объединить в т.н. \(6N\)-мерное фазовое пространство (ФП), каждая точка \(({\mathbf p},{\mathbf q})\) которого содержит все импульсы \({\mathbf p}=p_{x,1},p_{y,1},p_{z,1},\cdots,p_{x,N},p_{y,N},p_{z,N}\), и все координаты \({\mathbf q}=q_{x,1},q_{y,1},q_{z,1},\cdots,q_{x,N},q_{y,N},q_{z,N}\) всех частиц в один момент времени \(t\). Любая точка ФП полностью описывает текущее состояние системы, поэтому её можно интерпретировать как микросостояние, а ФП — как совокупность всех возможных микросостояний. По мере эволюции системы точки движутся в ФП, проделанные ими пути называют фазовыми траекториями. Плотности микросостояний в ФП приписывают функцию распределения \(\rho(\mathbf{p},\mathbf{q},t)\), она в состоянии равновесия не должна зависеть от времени, в этом случае точки распределены в ФП с неизменной плотностью; можно сформулировать и более общее утверждение, называемое теоремой Лиувилля: \[ \frac{d\rho(\mathbf{p},\mathbf{q},t)}{dt}=0, \] согласно которой плотность фазовых точек при их движении по фазовым траекториям постоянна. Если \(\rho(\mathbf{p},\mathbf{q})\) не зависит от времени явно, то \[ \frac{\partial \rho(\mathbf{p},\mathbf{q})}{\partial t}=0. \] Теорема Лиувилля является следствием уравнений Гамильтона: \[ \dot{q}_{i}=\frac{\partial H(\mathbf{p},\mathbf{q})}{\partial p_{i}} \quad \text{и} \quad \dot{p}_{i}=-\frac{\partial H(\mathbf{p},\mathbf{q})}{\partial q_{i}}, \] где \(H\) — гамильтониан системы, точкой обозначена полная производная по времени.

Если в результате МД было получено достаточно много наблюдений, можно сказать, что система прошла длинную фазовую траекторию и, вероятно, побывала во всех областях ФП, совокупность этих точек будет задавать геометрический прообраз функции распределения \(\rho(\mathbf{p},\mathbf{q})\). Отметим, что закон распределения \(\rho(\mathbf{p},\mathbf{q})\) полностью описывает ТД систему в статистическом смысле.

3 RMSD

Полученные траектории атомов можно анализировать с позиции конфигурации молекул. Важным параметром конкретной конфигурации (и всех точек траектории в целом) является RMSD² относительно эталонной структуры, по сути это среднеквадратичное отклонение координат атомов от координат соответствующих атомов в эталонной конформации: \[ \mathrm{RMSD}=\sqrt{\dfrac{1}{\Omega}\sum_{i}\omega_{i}\left(\mathbf{r}_{i}-\mathbf{r}_{i}^{\text{ref}}\right)^{2}}, \] где \(r_{i}\) и \(r_{i}^{\text{ref}}\) — координаты \(i\)го атома в текущей и эталонной конфигурации, \(\omega_{i}\) — вес этой пары (обычно принимается за 1), \(\Omega\) — сумма всех весов.

² Root-Mean-Square Deviation