diff --git a/ru/_posts/2017-12-25-practice-of-md-simulation-p1.md b/ru/_posts/2017-12-25-practice-of-md-simulation-p1.md new file mode 100644 index 0000000..f00923f --- /dev/null +++ b/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/), +что означает, что вы можете свободно распространять, модифицировать, создавать +производные произведения с условием, что вы укажите меня автором оригинального +произведения и, в том случае, если вы хотите распространять производную работу, +она должна быть под той же лицензией. + + + +## Определения метода молекулярной динамики + +Метод молекулярной динамики - метод, в котором временная эволюция системы +взаимодействующих атомов или частиц отслеживается интегрированием их уравнений +движения ([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\\) ансамбль). + +В качестве первого приближения можно считать, что объем можно хитрым образом +выразить через давление и температуру в системе (уравнение Менделеева-Клапейрона). +Тогда задача сводится к - так или иначе - заданию определенного давления и +температуры. Для поддержания определенной температуры в системе используется алгоритм, +который описывает термостат (обмен температуры системы с окружающей средой). Для +давления, соответственно, баростат - "обмен" давления с окружающей средой. + +## Молекулярно-динамические траектории + +Под траекторией мы будем подразумевать зависимость координат частиц от времени. +Также, в виду лени, сюда мы будем относить и все производные физические величины, +которые мы можем рассчитать напрямую в ходе моделирования - скорости, ускорения, +температура, давление, различные энергии и т.д.