Ограничение (вычислительная химия)

В вычислительной химии алгоритм ограничений — это метод удовлетворения ньютоновского движения твердого тела, состоящего из массовых точек. Алгоритм ограничения используется для обеспечения сохранения расстояния между массовыми точками. Общие шаги включают в себя: (i) выбрать новые неограниченные координаты (внутренние координаты), (ii) ввести явные ограничивающие силы, (iii) неявно минимизировать ограничивающие силы с помощью техники множителей Лагранжа или методов проекции.

Алгоритмы ограничений часто применяются для моделирования молекулярной динамики . Хотя такое моделирование иногда выполняется с использованием внутренних координат, которые автоматически удовлетворяют ограничениям по длине связи, валентному углу и углу кручения, моделирование также может выполняться с использованием явных или неявных ограничивающих сил для этих трех ограничений. Однако явные сдерживающие силы приводят к неэффективности; для получения траектории заданной длины требуется больше вычислительной мощности. Поэтому обычно предпочтительнее использовать внутренние координаты и средства решения неявных силовых ограничений.

Алгоритмы ограничений достигают вычислительной эффективности за счет пренебрежения движением по некоторым степеням свободы. Например, в атомистической молекулярной динамике обычно длина ковалентных связей с водородом ограничена; однако алгоритмы ограничений не следует использовать, если вибрации вдоль этих степеней свободы важны для изучаемого явления.

Математическая основа

Движение набора из N частиц можно описать системой обыкновенных дифференциальных уравнений второго порядка, вторым законом Ньютона, который можно записать в матричной форме

$\mathbf {M} \cdot {\frac {d^{2}\mathbf {q} }{dt^{2}}}=\mathbf {f} =-{\frac {\partial V}{\partial \mathbf {q} }}$

где M — матрица масс , а q — вектор обобщенных координат , описывающих положения частиц. Например, вектор q может представлять собой 3N декартовы координаты положений частиц r _k , где k изменяется от 1 до N ; в отсутствие ограничений M будет размером 3N x 3N диагональной квадратной матрицей масс частиц . Вектор f представляет обобщенные силы, а скаляр V ( q ) представляет потенциальную энергию, оба из которых являются функциями обобщенных координат q .

Если M присутствуют ограничений, координаты также должны удовлетворять M независимым от времени алгебраическим уравнениям.

$g_{j}(\mathbf {q} )=0$

индекс j изменяется от 1 до M. где Для краткости эти функции g _i сгруппированы в M -мерный вектор g ниже. Задача состоит в том, чтобы решить совокупность дифференциально-алгебраических уравнений (ДАУ), а не только обыкновенные дифференциальные уравнения (ОДУ) второго закона Ньютона.

Эту проблему детально изучил Жозеф Луи Лагранж , который изложил большинство методов ее решения. Самый простой подход — определить новые обобщенные координаты, которые не имеют ограничений; этот подход исключает алгебраические уравнения и снова сводит задачу к решению обыкновенного дифференциального уравнения. Такой подход используется, например, при описании движения твердого тела; положение и ориентация твердого тела могут быть описаны шестью независимыми неограниченными координатами, а не описывать положения составляющих его частиц и ограничения между ними, которые поддерживают их относительные расстояния. Недостатком этого подхода является то, что уравнения могут стать громоздкими и сложными; например, матрица масс M может стать недиагональной и зависеть от обобщенных координат.

Второй подход заключается во введении явных сил, которые поддерживают ограничение; например, можно ввести сильные пружинные силы, которые будут обеспечивать соблюдение расстояний между точками массы внутри «твердого» тела. Две трудности этого подхода заключаются в том, что ограничения не выполняются точно, а сильные взаимодействия могут потребовать очень коротких временных шагов, что делает моделирование неэффективным с вычислительной точки зрения.

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

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

Методы внутренних координат

Самый простой подход к удовлетворению ограничений в минимизации энергии и молекулярной динамике состоит в том, чтобы представить механическую систему в так называемых внутренних координатах, соответствующих неограниченным независимым степеням свободы системы. Например, двугранные углы белка представляют собой независимый набор координат, которые определяют положения всех атомов, не требуя каких-либо ограничений. Трудность таких подходов с использованием внутренних координат двояка: ньютоновские уравнения движения становятся намного более сложными, а внутренние координаты может быть трудно определить для циклических систем ограничений, например, при сморщивании колец или когда белок имеет дисульфидную связь.

Оригинальные методы эффективной рекурсивной минимизации энергии во внутренних координатах были разработаны Го и его коллегами.

Эффективные рекурсивные решатели ограничений по внутренней координате были распространены на молекулярную динамику. Аналогичные методы позже были применены и к другим системам.

Методы, основанные на множителях Лагранжа

Разрешение ограничений жесткой молекулы воды с использованием множителей Лагранжа : а) неограниченные положения получаются после определенного временного шага моделирования, б) вычисляются градиенты каждого ограничения по каждой частице и в) множители Лагранжа вычисляются для каждого градиента, такого как что ограничения удовлетворены.

В большинстве моделей молекулярной динамики, в которых используются алгоритмы ограничений, ограничения применяются с использованием метода множителей Лагранжа. Учитывая набор из n линейных ( голономных ) ограничений в момент времени t ,

$\sigma _{k}(t):=\|\mathbf {x} _{k\alpha }(t)-\mathbf {x} _{k\beta }(t)\|^{2}-d_{k}^{2}=0,\quad k=1\ldots n$

где $\scriptstyle \mathbf {x} _{k\alpha }(t)$ и $\scriptstyle \mathbf {x} _{k\beta }(t)$ — это положения двух частиц, участвующих в k -м ограничении в момент времени t и $d_{k}$ – предписанное расстояние между частицами.

Силы, возникающие из-за этих ограничений, добавляются в уравнения движения, в результате чего для каждой из N частиц в системе

${\frac {\partial ^{2}\mathbf {x} _{i}(t)}{\partial t^{2}}}m_{i}=-{\frac {\partial }{\partial \mathbf {x} _{i}}}\left[V(\mathbf {x} _{i}(t))-\sum _{k=1}^{n}\lambda _{k}\sigma _{k}(t)\right],\quad i=1\ldots N.$

Добавление сил ограничения не меняет общую энергию, поскольку чистая работа, совершаемая силами ограничения (взятая за набор частиц, на которые действуют ограничения), равна нулю. Обратите внимание, что знак на $\lambda _{k}$ произвольно и некоторые ссылки иметь противоположный знак.

Интегрируя обе части уравнения по времени, получаем ограниченные координаты частиц в данный момент времени: $t+\Delta t$, даны,

$\mathbf {x} _{i}(t+\Delta t)={\hat {\mathbf {x} }}_{i}(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1},\quad i=1\ldots N$

где ${\hat {\mathbf {x} }}_{i}(t+\Delta t)$ — это неограниченное (или нескорректированное) положение i -й частицы после интегрирования неограниченных уравнений движения.

Чтобы удовлетворить ограничения $\sigma _{k}(t+\Delta t)$ на следующем временном шаге множители Лагранжа должны быть определены как следующее уравнение:

$\sigma _{k}(t+\Delta t):=\left\|\mathbf {x} _{k\alpha }(t+\Delta t)-\mathbf {x} _{k\beta }(t+\Delta t)\right\|^{2}-d_{k}^{2}=0.$

Это подразумевает решение системы $n$ нелинейные уравнения

$\sigma _{j}(t+\Delta t):=\left\|{\hat {\mathbf {x} }}_{j\alpha }(t+\Delta t)-{\hat {\mathbf {x} }}_{j\beta }(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}\left(\Delta t\right)^{2}\left[{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{j\beta }}}m_{j\beta }^{-1}\right]\right\|^{2}-d_{j}^{2}=0,\quad j=1\ldots n$

одновременно для $n$ неизвестные множители Лагранжа $\lambda _{k}$.

Эта система $n$ нелинейные уравнения в $n$ неизвестные обычно решаются с использованием метода Ньютона – Рафсона , где вектор решения ${\underline {\lambda }}$ обновляется с помощью

${\underline {\lambda }}^{(l+1)}\leftarrow {\underline {\lambda }}^{(l)}-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}(t+\Delta t)$

где $\mathbf {J} _{\sigma }$ – якобиан уравнений σ _k :

$\mathbf {J} =\left({\begin{array}{cccc}{\frac {\partial \sigma _{1}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{1}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{1}}{\partial \lambda _{n}}}\\[5pt]{\frac {\partial \sigma _{2}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{2}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{2}}{\partial \lambda _{n}}}\\[5pt]\vdots &\vdots &\ddots &\vdots \\[5pt]{\frac {\partial \sigma _{n}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{n}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{n}}{\partial \lambda _{n}}}\end{array}}\right).$

Поскольку не все частицы вносят вклад во все ограничения, $\mathbf {J} _{\sigma }$ является блочной матрицей и может быть решена индивидуально до блочной единицы матрицы. Другими словами, $\mathbf {J} _{\sigma }$ можно решить индивидуально для каждой молекулы.

Вместо постоянного обновления вектора ${\underline {\lambda }}$, итерацию можно начать с ${\underline {\lambda }}^{(0)}=\mathbf {0}$, что приводит к более простым выражениям для $\sigma _{k}(t)$ и ${\frac {\partial \sigma _{k}(t)}{\partial \lambda _{j}}}$. В этом случае

$J_{ij}=\left.{\frac {\partial \sigma _{j}}{\partial \lambda _{i}}}\right|_{\mathbf {\lambda } =0}=2\left[{\hat {x}}_{j\alpha }-{\hat {x}}_{j\beta }\right]\left[{\frac {\partial \sigma _{i}}{\partial x_{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{i}}{\partial x_{j\beta }}}m_{j\beta }^{-1}\right].$

затем $\lambda$ обновляется до

$\mathbf {\lambda } _{j}=-\mathbf {J} ^{-1}\left[\left\|{\hat {\mathbf {x} }}_{j\alpha }(t+\Delta t)-{\hat {\mathbf {x} }}_{j\beta }(t+\Delta t)\right\|^{2}-d_{j}^{2}\right].$

После каждой итерации неограниченные положения частиц обновляются с использованием

${\hat {\mathbf {x} }}_{i}(t+\Delta t)\leftarrow {\hat {\mathbf {x} }}_{i}(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1}.$

Затем вектор сбрасывается на

${\underline {\lambda }}=\mathbf {0} .$

Описанная выше процедура повторяется до решения уравнений связей: $\sigma _{k}(t+\Delta t)$, сходится к заданному допуску на числовую ошибку.

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

Алгоритм SETTLE

Алгоритм SETTLE решает систему нелинейных уравнений аналитически относительно $n=3$ ограничения в постоянное время. Хотя он не масштабируется для большего количества ограничений, он очень часто используется для ограничения жестких молекул воды, которые присутствуют почти во всех биологических моделях и обычно моделируются с использованием трех ограничений (например, SPC/E и TIP3P модели воды ).

Алгоритм SHAKE

Алгоритм SHAKE был впервые разработан для удовлетворения ограничений геометрии связи во время моделирования молекулярной динамики. Затем метод был обобщен для обработки любых голономных ограничений, например тех, которые необходимы для поддержания постоянных валентных углов или молекулярной жесткости.

В алгоритме SHAKE система нелинейных уравнений ограничений решается с использованием метода Гаусса – Зейделя , который аппроксимирует решение линейной системы уравнений с использованием метода Ньютона – Рафсона ;

${\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}.$

Это равнозначно предположению, что $\mathbf {J} _{\sigma }$ является диагонально доминирующим и решает задачу $k$уравнение только для $k$ неизвестный. На практике мы вычисляем

${\begin{aligned}\lambda _{k}&\leftarrow {\frac {\sigma _{k}(t)}{\partial \sigma _{k}(t)/\partial \lambda _{k}}},\\[5pt]\mathbf {x} _{k\alpha }&\leftarrow \mathbf {x} _{k\alpha }+\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{k\alpha }}},\\[5pt]\mathbf {x} _{k\beta }&\leftarrow \mathbf {x} _{k\beta }+\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{k\beta }}},\end{aligned}}$

для всех $k=1\ldots n$ итеративно до тех пор, пока уравнения ограничений $\sigma _{k}(t+\Delta t)$ решаются с заданным допуском.

Стоимость расчета каждой итерации составляет ${\mathcal {O}}(n)$, а сами итерации сходятся линейно.

Неитеративная форма SHAKE была разработана позже.

Существует несколько вариантов алгоритма SHAKE. Хотя они различаются в том, как они вычисляют или применяют сами ограничения, ограничения по-прежнему моделируются с использованием множителей Лагранжа, которые вычисляются с использованием метода Гаусса – Зейделя.

Оригинальный алгоритм SHAKE способен ограничивать как жесткие, так и гибкие молекулы (например, воду, бензол и бифенил ) и вносит незначительную ошибку или дрейф энергии в моделирование молекулярной динамики. Одна из проблем SHAKE заключается в том, что количество итераций, необходимых для достижения определенного уровня сходимости, действительно увеличивается по мере усложнения молекулярной геометрии. Чтобы достичь 64-битной компьютерной точности (относительная допуск $\approx 10^{-16}$) в типичном молекулярно-динамическом моделировании при температуре 310K модель воды с 3 участками, имеющая 3 ограничения для поддержания молекулярной геометрии, требует в среднем 9 итераций (что составляет 3 на участок за временной шаг). с 12 узлами Модель бутана с 4 узлами и 5 ограничениями требует 17 итераций (22 на сайт), модель бензола с 6 узлами с 12 ограничениями требует 36 итераций (72 на сайт), а модель бифенила с 29 ограничениями требует 92 итераций ( 229 на сайт за временной шаг). Следовательно, требования к процессору алгоритма SHAKE могут стать значительными, особенно если молекулярная модель имеет высокую степень жесткости.

Более позднее расширение метода QSHAKE ( Quaternion SHAKE) было разработано как более быстрая альтернатива для молекул, состоящих из жестких единиц, но оно не является универсальным. Он удовлетворительно работает для жестких петель, таких как системы ароматических колец , но QSHAKE не работает для гибких петель, например, когда белок имеет дисульфидную связь.

Дальнейшие расширения включают RATTLE, ВИГГЛ, и МШАКЕ.

Хотя RATTLE работает так же, как SHAKE, однако, используя схему интегрирования времени Velocity Verlet , WIGGLE расширяет SHAKE и RATTLE, используя начальную оценку множителей Лагранжа. $\lambda _{k}$ в зависимости от скорости частиц. Стоит отметить, что MSHAKE вычисляет поправки к силам ограничения , обеспечивая лучшую сходимость.

Последней модификацией алгоритма SHAKE является алгоритм P-SHAKE. это применяется к очень жестким или полужестким молекулам . P-SHAKE вычисляет и обновляет предварительный кондиционер, который применяется к градиентам ограничений перед итерацией SHAKE, вызывая Якобиан $\mathbf {J} _{\sigma }$ стать диагональным или сильно диагонально доминирующим. Таким образом, разделенные ограничения сходятся гораздо быстрее (квадратично, а не линейно) за счет ${\mathcal {O}}(n^{2})$.

Алгоритм M-SHAKE

Алгоритм M-SHAKE решает нелинейную систему уравнений напрямую методом Ньютона . На каждой итерации линейная система уравнений

${\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}$

решается точно с использованием LU-разложения . Каждая итерация стоит ${\mathcal {O}}(n^{3})$ операций, однако решение сходится квадратично , требуя меньше итераций, чем SHAKE.

Это решение было впервые предложено в 1986 году Чиккотти и Рикартом. под названием «матричный метод», однако отличался решением линейной системы уравнений. Чиккотти и Рикарт предлагают инвертировать матрицу $\mathbf {J} _{\sigma }$ напрямую, но делаем это только один раз, в первой итерации. Тогда первая итерация стоит ${\mathcal {O}}(n^{3})$ операций, тогда как последующие итерации стоят всего ${\mathcal {O}}(n^{2})$ операции (для умножения матрицы на вектор). Однако за это улучшение приходится платить, поскольку якобиан больше не обновляется, сходимость является только линейной , хотя и с гораздо большей скоростью, чем для алгоритма SHAKE.

Несколько вариантов этого подхода, основанного на методах разреженной матрицы, были изучены Barth et al. .

Алгоритм ФОРМА

Алгоритм ФОРМА является многоцентровым аналогом SHAKE для ограничения твердых тел тремя и более центрами. Как и в случае с SHAKE, делается неограниченный шаг, а затем он корректируется путем непосредственного расчета и применения матрицы вращения твердого тела, которая удовлетворяет:

$L^{\text{rigid}}\left(t+{\frac {\Delta t}{2}}\right)=L^{\text{nonrigid}}\left(t+{\frac {\Delta t}{2}}\right)$

Этот подход включает в себя одну диагонализацию матрицы 3×3, за которой следуют три или четыре быстрых итерации Ньютона для определения матрицы вращения. SHAPE обеспечивает ту же траекторию, что и полностью конвергентная итеративная SHAKE, однако она оказывается более эффективной и точной, чем SHAKE, когда применяется к системам, включающим три или более центров. Он расширяет возможности ограничений, подобных SHAKE, на линейные системы с тремя и более атомами, плоские системы с четырьмя и более атомами, а также на значительно более крупные жесткие структуры, где SHAKE неразрешим. Он также позволяет связывать твердые тела с одним или двумя общими центрами (например, пептидными плоскостями) путем итеративного решения ограничений твердого тела тем же основным способом, которым SHAKE используется для атомов, включающих более одного ограничения SHAKE.

Алгоритм LINCS

Альтернативный метод ограничений, LINCS (линейный решатель ограничений), был разработан в 1997 году Хессом, Беккером, Берендсеном и Фрайе. и был основан на методе Эдберга, Эванса и Морриса 1986 года (EEM), и его модификация Бараньяи и Эванса (БЕЛАРУСЬ).

LINCS применяет множители Лагранжа к силам ограничения и находит множители, используя разложение в ряд для аппроксимации обратного якобиана. $\mathbf {J} _{\sigma }$:

$(\mathbf {I} -\mathbf {J} _{\sigma })^{-1}=\mathbf {I} +\mathbf {J} _{\sigma }+\mathbf {J} _{\sigma }^{2}+\mathbf {J} _{\sigma }^{3}+\cdots$

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

Сообщается, что LINCS работает в 3–4 раза быстрее, чем SHAKE.

Гибридные методы

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

См. также

• Молекулярная динамика
• Программное обеспечение для моделирования молекулярной механики