add md1 first edition

ru/_posts/2017-12-25-practice-of-md-simulation-p1.md

@@ -0,0 +1,253 @@
+---
+category: ru
+type: paper
+hastr: false
+layout: paper
+tags: molecular dynamics, gromacs
+title: Практическая молекулярная динамика. Часть 1
+short: practice-of-md-simulation-p1
+use_math: true
+---
+Серия статей, в которых я постараюсь научить людей работать с методом классической
+молекулярной динамики и программным пакетом GROMACS в частности. Они основаны на
+моих приватных лекциях студентам, моем собственном опыте и обсуждениях с коллегами.
+Я постараюсь рассказывать в доступной форме специально для неподготовленных людей,
+таким образом (я надеюсь), может рассматриваться как обучающее руководство для
+новичков. Любые вопросы, комментарии и замечания приветствуются.
+
+Часть 1 содержит введение и некоторые основы метода молекулярной динамики.
+
+Все статьи лицензированы под [CC BY-SA](https://creativecommons.org/licenses/by-sa/4.0/),
+что означает, что вы можете свободно распространять, модифицировать, создавать
+производные произведения с условием, что вы укажите меня автором оригинального
+произведения и, в том случае, если вы хотите распространять производную работу,
+она должна быть под той же лицензией.
+
+<!--more-->
+
+## Определения метода молекулярной динамики
+
+Метод молекулярной динамики - метод, в котором временная эволюция системы
+взаимодействующих атомов или частиц отслеживается интегрированием их уравнений
+движения ([Wikipedia](https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%BA%D0%BB%D0%B0%D1%81%D1%81%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%BE%D0%B9_%D0%BC%D0%BE%D0%BB%D0%B5%D0%BA%D1%83%D0%BB%D1%8F%D1%80%D0%BD%D0%BE%D0%B9_%D0%B4%D0%B8%D0%BD%D0%B0%D0%BC%D0%B8%D0%BA%D0%B8)). На
+практике это означает примерно следующее
+
+* Для описания движения используются положения классической ньютоновской механики
+* Силы взаимодействия между частицами задаются в форме каких либо классических
+потенциальных функций.
+
+Давайте рассмотрим эти положения подробнее.
+
+Надеюсь, вы знаете что такое ньютоновская механика, это что то, строящееся на
+втором законе Ньютона:
+
+$$
+m\ddot{x} = F
+$$
+
+идея данного положения заключается в том, что у нас есть некоторая сила, которая
+задается по потенциалу (ака первая производная от энергии), есть известная масса,
+на основании которых мы можем рассчитать ускорение. В дальнейшем простым
+(на самом деле нет) интегрированием можно получить скорости и - в дальнейшем -
+координаты.
+
+Возможно, для понимания было бы проще использовать уравнение Лагранжа, которое
+выглядит следующим образом:
+
+$$
+\frac{d}{dt} \frac{\partial L}{\partial \dot{x_i}} - \frac{\partial L}{\partial x_i} = Q^n_i
+$$
+
+(ну или другие варианты, которые будут рассмотрены ниже)
+
+Второе положение по счастливой случайности описывает то, каким образом мы можем
+получить значение силы, действующий на частицу (атом) в заданный момент времени.
+В рамках молекулярной механики используется в основном четыре типа взаимодействия:
+
+* внутримолекулярные (они же растяжение связи и деформация _плоских_ углов),
+  которые описываются гармоническим потенциалом:
+
+  $$  
+  U = \frac{k}{2}(r-r_0)^2
+  $$
+
+* внутримолекулярные для описания торсионных углов (там есть свои нюансы, но
+  описываются они, к сожалению, уже не гармониками);
+* межмолекулярные, описываемые кулоновским законом;
+* межмолекулярные (они же поляризационные, они же ван-дер-ваальсовые).
+
+Таким образом, энергия взаимодействия одной частицы со всей системой описывается
+следующим уравнением:
+
+$$
+U_{total} = U_{12} + U_{13} + U_{14} + U_{vdw} + U_{coul}
+$$
+
+где \\(U_{12}\\) и \\(U_{13}\\) - гармоники, описывающие взаимодействие между
+соседними частицами (оно же растяжение связей) и через одну частицу (оно же деформация
+плоских углов) соответственно, \\(U_{14}\\) - деформация торсионных углов,
+\\(U_{vdw}\\) и \\(U_{coul}\\) - энергии ван-дер-ваальсового и кулоновского
+взаимодействия соответственно.
+
+Если обобщить, то потенциальная энергия системы описывается некоторой функцией,
+которую можно однозначно определить имея текущие координаты частицы и координаты
+других частиц в системе:
+
+$$
+U_i = f(r_i, \\{r_j\\}^N_0)
+$$
+
+Ну а полная энергия системы будет суммой потенциальных и кинетических энергий
+всех частиц в системе, что-то вроде:
+
+$$
+E = \sum_{i}{U_i} + \sum_{i}{\frac{m_i {\dot x}_i^2}{2}}
+$$
+
+То есть, зная координаты системы в заданный момент времени,
+можно рассчитать силу и энергию, на основании которых можно некоторым магическим
+образом рассчитать новые координаты.
+
+## Интегрирование уравнений движения
+
+По сути, интегрирование уравнений движения есть решение уравнений Ньютона (или
+Лагранжа, в зависимости от того, в рамках какого формализма мы рассматриваем
+систему). Поскольку машины глупые, а вычислительные мощности ограничены, то для
+их решения используются приближенные методы.
+
+### Методы численного интегрирования
+
+Для численного (в отличие от аналитического с первообразными функциями, знакомого
+по школьной программе) интегрирования используется геометрическое определение
+интегрирования - определенный интеграл от заданной функции есть площадь фигуры под
+ней.
+
+Для того, чтобы рассчитать площадь фигуры, можно разбить ее на маленькие части,
+например, прямоугольные трапеции с заданной шириной \\(\Delta x\\). Площадь такой
+трапеции можно рассчитать по формуле
+
+$$
+S_i = \Delta x \frac{f(x_i) + f(x_i + \Delta x)}{2}
+$$
+
+тогда площадь всей фигуры, если ее разбить на \\(N\\) равных по ширине трапеций,
+будет равна
+
+$$
+\int f(x) dx = \sum_{i=1}^{N} {S_i} = \sum_{i=0}^{N-1} { \Delta x \frac{f(x + i \Delta x) + f(x + (i + 1) \Delta x)}{2} }
+$$
+
+В численных методах (и дальше по тексту) \\(\Delta x\\) мы будем называть шагом
+интегрирования и будем считать его предельно малой величиной.
+
+Следует понимать, что приведенный здесь метод является, пожалуй, самым простым
+для понимания, но отнюдь не единственным (и не самым точным). Более подробно
+вопросы численного интегрирования будут рассмотрены дальше.
+
+### Алгоритм Верле
+
+Допустим, у нас задано время \\(t\\) и малый шаг интегрирования \\(\Delta t\\).
+Тогда в моменты времени \\(t + \Delta t\\) и \\(t - \Delta t\\) координаты можно
+задать как нечто, что раскладывается в ряд Тейлора в окрестности \\(t\\)
+(пренебрегая членами порядка 3 и выше, поскольку мы предполагаем, что шаг
+интегрирования бесконечно мал):
+
+$$
+x(t + \Delta t) = x(t) + \dot x (t) \Delta t + \frac{\ddot x(t) {\Delta t}^2}{2}
+$$
+
+$$
+x(t - \Delta t) = x(t) - \dot x (t) \Delta t + \frac{\ddot x(t) {\Delta t}^2}{2}
+$$
+
+Ничто нам не мешает сложить два уравнения друг с другом и получить следующее:
+
+$$
+x(t + \Delta t) + x(t - \Delta t) = 2x(t) + \ddot x(t) {\Delta t}^2
+$$
+
+или
+
+$$
+x(t + \Delta t) = 2x(t) - x(t - \Delta t) + \ddot x(t) {\Delta t}^2
+$$
+
+Если перевести на человеческий язык. то, зная координаты частицы в момент времени
+\\(t\\) и немного (на \\(\Delta t\\)) раньше, а также зная ускорение, мы можем
+рассчитать координаты частицы в будущем на \\(\Delta t\\).
+
+Ускорение можно рассчитать исходя из силы, используя второй закон Ньютона (которую,
+в свою очередь, легко рассчитать дифференцированием энергии), координаты частиц
+мы, вроде бы, знаем. Но есть нюанс. Безусловно, мы знаем координаты частиц в
+текущий момент времени, но откуда взять координаты частиц из прошлого, если
+текущий момент времени равен 0?
+
+Для определения координат частиц в предыдущий 0 момент времени используется
+примерно следующий хак. Допустим, мы проводим моделирование при заданной температуре
+\\(T\\); воспользовавшись аппаратом статистической физики, мы можем описать
+распределение частиц по скоростям (например,
+[Максвелловское распределение](https://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D0%9C%D0%B0%D0%BA%D1%81%D0%B2%D0%B5%D0%BB%D0%BB%D0%B0)). Таким образом, можно случайным образом задать частице такую
+скорость, чтобы температура всей системы была равна заданной. Далее, рассчитать
+ускорение частицы в нулевой момент времени и, в дальнейшем, собственно координаты
+частиц в момент времени \\(t - \Delta t\\).
+
+inb4: Следует понимать, что указанный алгоритм имеет свои недостатки. Например,
+поскольку скорости явно не учитываются при интегрировании уравнений движения,
+их расчет имеет некоторую погрешность, что приводит к значительно амплитуде
+изменения скорости одной частицы со временем (и, как следствие, погрешность при
+определении кинетической энергии). Также, указанный алгоритм довольно упрощен,
+для уточнения конкретных реализаций следует читать исходный код (если доступен)
+соответствующих программных пакетов или научные статьи.
+
+## Силовые поля
+
+Внимательный (лол) читатель мог обратить внимание, что, даже если рассматривать
+относительно простые случаи потенциальных энергий - энергии гармонических
+осцилляторов - они содержат некоторые магические параметры - \\(k\\) и \\(r_0\\).
+Если их физический смысл - я надеюсь - не вызывает вопросов - \\(k\\) - это что-то
+похожее на жесткость пружины, а \\(r_0\\) - расстояние в точке минимума - то вопрос,
+откуда их брать, может вызвать некоторые сложности.
+
+Для их определения, обычно решается обратная задача - исходя из известной энергии
+рассчитываются вторые производные (они же \\(k\\)), а \\(r_0\\) рассчитывается
+исходя из минимума потенциальной энергии. В качестве опорных значений для расчетов
+энергии используются _неэмерические_ методы (в противовес _эмпирической_ классической
+молекулярной механике), поголовно состоящие из квантово-химических методов.
+
+Совокупность же параметров, рассчитанных в рамках одного метода с усреднением
+параметров по определенному принципу (равно как и подгон значений под определенные
+вполне материальные физические величины) в дальнейшем мы будем называть силовым
+полем. Вопрос выбора того или иного силового поля, а также принципы, которые лежат
+в основе их определения, мы рассмотрим позже.
+
+## Термостат, баростат и все-все-все
+
+В методах статистической физики используются вполне определенные _ансабмли_. У
+них есть хитрые названия, которые я использовать не буду, потому что они не очень
+понятны. Принцип их образования примерно следующий: берутся определенные физические
+величины, которые нам известны в заданной конфигурации, а остальные физические
+величины рассчитываются исходя из определенных.
+
+В рамках данных статей я буду рассматривать только те системы, в которых число
+частиц постоянно. Некоторые величины, например, энергия, довольно трудно
+зафиксировать, поэтому они также не используются при моделировании. По факту,
+как правило, фиксируется две из трех величин: давление, температура и объем:
+
+* канонический ансамбль, в котором фиксируется число частиц, объем и температура
+(он же \\(NVT\\) ансамбль);
+* изотермо-изобарический ансамбль, в котором фиксируется число частиц, температура
+и давление (он же \\(NpT\\) ансамбль).
+
+В качестве первого приближения можно считать, что объем можно хитрым образом
+выразить через давление и температуру в системе (уравнение Менделеева-Клапейрона).
+Тогда задача сводится к - так или иначе - заданию определенного давления и
+температуры. Для поддержания определенной температуры в системе используется алгоритм,
+который описывает термостат (обмен температуры системы с окружающей средой). Для
+давления, соответственно, баростат - "обмен" давления с окружающей средой.
+
+## Молекулярно-динамические траектории
+
+Под траекторией мы будем подразумевать зависимость координат частиц от времени.
+Также, в виду лени, сюда мы будем относить и все производные физические величины,
+которые мы можем рассчитать напрямую в ходе моделирования - скорости, ускорения,
+температура, давление, различные энергии и т.д.