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

Алгоритм Верле — это численный метод интегрирования уравнений движения, используемый в молекулярной динамике и других областях вычислительной физики для расчёта траекторий частиц во времени. Алгоритм относится к классу симплектических интеграторов, что означает сохранение фазового объёма и энергии системы в долгосрочной перспективе при отсутствии диссипативных сил. Разработан в 1960-х годах французским физиком Лу Верле (Loup Verlet) для моделирования динамики атомов и молекул.

История

Метод был предложен Лу Верле в 1967 году в статье «Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules» (опубликована в журнале Physical Review). Верле использовал алгоритм для моделирования жидкого аргона с потенциалом Леннарда-Джонса, что стало одним из первых успешных применений компьютерного моделирования в статистической физике. До появления алгоритма Верле в молекулярной динамике применялись менее точные методы, такие как алгоритм Эйлера, которые приводили к быстрому дрейфу энергии.

В 1980-х годах алгоритм был усовершенствован в виде «скоростной версии» (velocity Verlet), предложенной Уильямом Свопом (William Swope) и коллегами, что позволило одновременно вычислять координаты и скорости частиц без потери симплектических свойств.

Математическая формулировка

Алгоритм Верле основан на разложении траектории частицы в ряд Тейлора до второго порядка. Для системы частиц с массами \( m_i \), координатами \( \mathbf{r}_i(t) \) и ускорениями \( \mathbf{a}_i(t) = \frac{\mathbf{F}_i(t)}{m_i} \) (где \( \mathbf{F}_i \) — сила, действующая на частицу) стандартная форма алгоритма записывается как:

\[ \mathbf{r}_i(t + \Delta t) = 2\mathbf{r}_i(t) - \mathbf{r}_i(t - \Delta t) + \mathbf{a}_i(t) (\Delta t)^2 + O(\Delta t^4) \]

где \( \Delta t \) — шаг интегрирования. Ошибка на каждом шаге составляет \( O(\Delta t^4) \) для координат, что делает алгоритм точным для малых \( \Delta t \). Скорости вычисляются через центральную разность:

\[ \mathbf{v}_i(t) = \frac{\mathbf{r}_i(t + \Delta t) - \mathbf{r}_i(t - \Delta t)}{2\Delta t} + O(\Delta t^2) \]

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

Скоростная версия (Velocity Verlet)

Скоростной алгоритм Верле (Velocity Verlet) решает эту проблему, вычисляя координаты и скорости на одном шаге:

1. Обновление координат:

\[ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \mathbf{v}_i(t) \Delta t + \frac{1}{2} \mathbf{a}_i(t) (\Delta t)^2 \]

1. Вычисление новых сил \( \mathbf{F}_i(t + \Delta t) \) на основе обновлённых координат.

1. Обновление скоростей:

\[ \mathbf{v}_i(t + \Delta t) = \mathbf{v}_i(t) + \frac{1}{2} \left[ \mathbf{a}_i(t) + \mathbf{a}_i(t + \Delta t) \right] \Delta t \]

Этот метод сохраняет симплектичность и имеет ошибку \( O(\Delta t^2) \) для скоростей, но при этом ошибка для координат остаётся \( O(\Delta t^4) \). Скоростная версия широко используется в современных пакетах молекулярной динамики, таких как GROMACS, LAMMPS, NAMD.

Свойства

Симплектичность

Алгоритм Верле является симплектическим интегратором, то есть сохраняет каноническую структуру фазового пространства (симплектическую форму). Это свойство обеспечивает долговременную стабильность энергии системы: при отсутствии внешних диссипативных сил полная энергия колеблется около постоянного значения без дрейфа. В отличие от несимплектических методов (например, Рунге — Кутты 4-го порядка), алгоритм Верле не приводит к экспоненциальному накоплению ошибок энергии.

Обратимость во времени

Алгоритм обратим во времени: если изменить знак шага \( \Delta t \), система вернётся в исходное состояние с точностью до ошибок округления. Это свойство важно для моделирования обратимых процессов, таких как классическая динамика в консервативных системах.

Устойчивость

Алгоритм устойчив для шагов \( \Delta t \), меньших характерного периода колебаний системы (например, для молекулярной динамики шаг обычно составляет 1–2 фемтосекунды). При превышении критического шага (порядка периода наивысшей частоты) возникает численная неустойчивость, приводящая к расходимости траекторий.

Применение

Молекулярная динамика

Алгоритм Верле является стандартным методом интегрирования в молекулярной динамике (МД). Он используется для моделирования движения атомов и молекул в жидкостях, твёрдых телах, биополимерах (белках, ДНК) и наноматериалах. В МД шаг интегрирования выбирается из условия разрешения самых быстрых колебаний (например, колебаний связей C–H, период которых около 10 фемтосекунд). Для увеличения шага применяются методы ограничения связей (например, алгоритм SHAKE, RATTLE), которые совместимы с алгоритмом Верле.

Астрофизика и небесная механика

В симуляциях гравитационных систем (например, моделирование движения планет, звёздных скоплений, галактик) алгоритм Верле используется для интегрирования уравнений Ньютона. Благодаря симплектичности он хорошо подходит для долгосрочных расчётов (миллионы шагов), где важно сохранение энергии и момента импульса.

Физика конденсированного состояния

Метод применяется для моделирования свойств материалов: фазовых переходов, диффузии, теплопроводности, механических деформаций. Например, в моделировании кристаллов с дефектами или аморфных структур.

Биофизика и фармакология

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

Ограничения

Необходимость хранения предыдущих координат

Стандартная версия алгоритма требует хранения координат на предыдущем шаге \( \mathbf{r}_i(t - \Delta t) \), что увеличивает объём памяти. Скоростная версия этого недостатка лишена, так как использует только текущие координаты и скорости.

Ошибка для скоростей

В стандартной версии скорости вычисляются с точностью \( O(\Delta t^2) \), что может быть недостаточно для задач, требующих высокой точности кинетической энергии (например, при расчёте температурных флуктуаций). В таких случаях предпочтительна скоростная версия.

Чувствительность к шагу

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

Варианты и модификации

Алгоритм «прыжок лягушки» (Leapfrog)

Эквивалентная форма алгоритма Верле, в которой скорости вычисляются на полуцелых шагах: \[ \mathbf{v}_i(t + \Delta t/2) = \mathbf{v}_i(t - \Delta t/2) + \mathbf{a}_i(t) \Delta t \] \[ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \mathbf{v}_i(t + \Delta t/2) \Delta t \] Этот вариант обладает теми же симплектическими свойствами, но требует меньше операций памяти.

Алгоритм Верле с термостатом

Для моделирования в каноническом ансамбле (NVT) алгоритм комбинируется с термостатами, такими как термостат Нозе — Гувера или термостат Ланжевена. В этом случае в уравнения добавляются стохастические или диссипативные члены, что нарушает симплектичность, но сохраняет статистическую точность.

Алгоритм Верле для уравнений с трением

В задачах с вязким трением (например, моделирование броуновского движения) используется модификация, включающая демпфирующий член. Пример — алгоритм Верле с вязкостью, применяемый в методе диссипативной динамики частиц (DPD).

Интересные факты

• Алгоритм Верле был независимо открыт несколькими исследователями: в 1967 году Лу Верле, а также ранее в 1960-х годах американским физиком Джоном Банчем (John Bunch) для решения задач небесной механики.
• В 1982 году Уильям Своп и коллеги опубликовали работу, в которой показали, что скоростная версия алгоритма является симплектической и имеет лучшие свойства по сравнению с исходной версией.
• Алгоритм Верле используется в симуляциях, охватывающих до миллиарда частиц (например, в моделировании коллоидных систем или в астрофизических N-телах).
• В 2013 году алгоритм был реализован на графических процессорах (GPU) для ускорения расчётов в молекулярной динамике, что позволило моделировать системы с миллионами атомов в реальном времени.

Источники

• Verlet, L. (1967). «Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules». Physical Review, 159(1), 98–103.
• Swope, W. C., Andersen, H. C., Berens, P. H., & Wilson, K. R. (1982). «A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters». The Journal of Chemical Physics, 76(1), 637–649.
• Frenkel, D., & Smit, B. (2002). «Understanding Molecular Simulation: From Algorithms to Applications». Academic Press.
• Allen, M. P., & Tildesley, D. J. (2017). «Computer Simulation of Liquids». Oxford University Press.
• Leimkuhler, B., & Reich, S. (2004). «Simulating Hamiltonian Dynamics». Cambridge University Press.