Молекулярная динамика (МД) - чрезвычайно популярный метод, который, помимо прочего, используется для моделирования движения атомов или молекул. МД приняли во многих областях, включая микромеханику, биофизику и химию. Одна из причин междисциплинарного характера MD - его относительная простота и надежность. Существует большое количество библиотек с открытым исходным кодом или свободного доступа для MD, включая LAMMPS, MDAnalysis и MDTraj, что делает проектирование и реализацию MD столь же простыми, как монтаж библиотеки и описание параметров моделирования.

Несмотря на легкодоступные библиотеки с открытым исходным кодом, чрезвычайно важно понимать истоки MD как метода моделирования, чтобы понять его результаты, построить модели и понять ограничения MD. Не вдаваясь в подробности численного вывода MD (есть много отличных учебников), я помогу вывести и продемонстрировать один из самых популярных алгоритмов интеграции положения: Velocity Verlet.

Что такое «интеграция позиций»?

Проще говоря, интегрирование положения в MD - это процесс, используемый для определения направления и расстояния, на которое движутся частицы с заданным временным шагом. Знание того, как группы частиц / молекул / положения атомов меняются с течением времени при наличии некоторого внешнего воздействия (например, силы, электрического поля), обычно может предоставить всю необходимую физическую информацию для понимания структурных, термодинамических или текучих свойств.

Истоки интеграции Velocity Verlet

Нашей отправной точкой будет обобщенный график, показывающий положение частицы (r) относительно времени (t).

Что нас интересует, так это то, как экстраполировать следующее положение частицы, учитывая всю информацию об остальной части системы (например, внешний стимул, положения всех других частиц и т. Д.). Короче говоря, нам нужно найти связь между r (t) и r (t + dt), где dt - временной шаг между итераций. Если мы предположим, что шаг dt между t и t + dt мал, мы можем аппроксимировать выражение для r (t + dt), используя расширение Тейлора:

O (dt⁴) - это термин значение ошибки, которое включает в себя завершающую ошибку из-за игнорирования dt⁴, dt⁵, dt⁶ и др. вклады. Здесь «над точками» обозначают первую, вторую и третью производные по времени по положению. Более привычными выражениями для этой записи могут быть следующие:

Для простоты я буду придерживаться обозначения над точкой во время вывода Velocity Verlet до окончательного выражения.

Возвращаясь к формуле. 1, мы можем изменить уравнение, чтобы выразить конечную разницу в положении между r (t) - ›r (t + dt):

Это также называется выражением «Вперед Эйлера», поскольку производная положения вычисляется при перемещении от r (t) - ›r (t + dt). И наоборот, мы также можем использовать выражение «Обратное Эйлера», вычислив производную, двигаясь от r (t-dt) - ›r (t).

Затем, если мы вычтем уравнение. 3 из уравнения. 2 (Уравнение 3 - Уравнение 2):

а затем упростим, решив для r (t + dt), мы получим:

Уравнение 4 - выражение для интегратора положения Верле, который является алгоритмом обновления положения. Однако это не самозапуск. Если вы представили симуляцию, начинающуюся с t = 0, для начала симуляции необходимо знать позицию в t = -dt.

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

Уравнение 5 можно переставить так, чтобы найти r (t-dt), которое принимает форму:

Затем мы можем вставить указанное выше для r (t-dt) в интегратор положения Верле (уравнение 4) и упростить его до формы:

Вышеприведенное выражение является окончательной формой интегратора положения Velocity Verlet. Заменив обозначение над точкой на более знакомые термины, мы получим:

Здесь мы видим, что все, что нам нужно знать для обновления положения, - это текущее положение, скорость, сила, масса и временной шаг. В отличие от интегратора Верле, Velocity Verlet запускается автоматически, потому что для его продвижения вперед не требуется никаких предварительных знаний о предыдущем состоянии системы.

Интегратор скорости для алгоритма Верле-Верле можно извлечь, сначала подставив (t + dt) вместо (t) членов в уравнении. 5 и уравнение 4, а затем подставив в уравнение. 4 для члена r (t + 2dt) в уравнении. 5, у нас осталось:

Или эквивалентная форма, которая часто используется:

Уравнения 6 и 7 - это интеграторы Верле скорости для положения и скорости. Как вы могли заметить, позиция r (t + dt) должна быть определена до определения скорости v (t + dt). Тогда только разумно, что обновление позиции должно произойти до обновления скорости. Имея это в виду, существует традиционная последовательность шагов, которые обычно используются при выполнении моделирования молекулярной динамики. Что-то вроде:

  1. Обновить все позиции
  2. Обновить все силы
  3. Обновить все скорости
  4. Применение граничных условий (включая контроль температуры и давления)
  5. Переместить глобальное время вперед на временной шаг (dt)
  6. Рассчитайте любой результат
  7. Повторите шаги 1–6 столько раз, сколько необходимо.

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

Мы можем продемонстрировать простое одномерное применение алгоритма Velocity Verlet, используя систему двух частиц, которые взаимодействуют с пружинным потенциалом (код представлен ниже)

Здесь пружина имеет расстояние покоя 5, но массы изначально разделены расстоянием в 4. Если вы присмотритесь, вы увидите, что максимальное расширение, кажется, увеличивается с каждым циклом ... так что есть ошибка в реализация / кодирование алгоритма интеграции…? Не обязательно…

Ошибка дискретизации

Помните все те термины об ошибках, которые игнорировались? (см. уравнение 1–5) На каждом этапе интегрирования накапливается небольшая ошибка, и со временем эта ошибка накапливается. Этот тип ошибки из-за усечения во время процесса интеграции называется ошибкой дискретизации. При использовании интегрирования Верле по скорости ошибка дискретизации обычно является экспоненциальной функцией временного шага (dt), что означает, что меньший временной шаг уменьшает общую накопленную ошибку. Это один из первых из многих компромиссов между точностью и вычислительной эффективностью при моделировании в MD. Использование чрезвычайно малого временного шага может привести к точным результатам, но за счет более длительного времени вычислений или сокращения общего времени моделирования, что может снизить статистическую надежность. Более крупный временной шаг может помочь в больших или более длительных временных масштабах, но за счет физической точности.

Если вам все еще не нравится зависимость накопления ошибок от шага по времени, эффект можно показать, запустив такое же моделирование, как указано выше, но с немного большими и меньшими шагами по времени. Поскольку исходное моделирование было dt = 0,01, я буду использовать dt = 0,02 и dt = 0,005 для большего и меньшего временных шагов.

Без какого-либо численного анализа совершенно очевидно, что выбранный временной шаг может оказать значительное влияние на результаты моделирования! Альтернативный способ продемонстрировать то, что мы можем визуально наблюдать на гифках выше, - это построить график эволюции полной энергии системы во времени. Поскольку мы имеем дело с пружиной в замкнутой среде, нам нужно, чтобы общая энергия (кинетическая энергия + потенциал) сохранялась с течением времени. В действительности мы, вероятно, увидим что-то вроде иллюстрации ниже.

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

Код для демонстрации MD Spring и сохранения gif

Спасибо за чтение - любые предложения и исправления приветствуются!

использованная литература

  1. Френкель, Даан и Беренд Смит. Понимание молекулярного моделирования от алгоритмов до приложений. Академик Пресс, 2002.
  2. Браун, Ефрем и др. «Лучшие практики для основ в молекулярном моделировании [статья v1.0]». Living Journal of Computational Molecular Science, vol. 1, вып. 1, 2019, DOI: 10.33011 / livecoms.1.1.5957.