25 дек 2017
Тeги: moleculardynamics,gromacs
Похожие посты:
Как загрузить скриншот в S3 с помощью linux
Awesome Widgets - Произвольные форматеры и макросы
Серия статей, в которых я постараюсь научить людей работать с методом классической молекулярной динамики и программным пакетом GROMACS в частности. Они основаны на моих приватных лекциях студентам, моем собственном опыте и обсуждениях с коллегами. Я постараюсь рассказывать в доступной форме специально для неподготовленных людей, таким образом (я надеюсь), может рассматриваться как обучающее руководство для новичков. Любые вопросы, комментарии и замечания приветствуются.
Часть 1 содержит введение и некоторые основы метода молекулярной динамики.
Все статьи лицензированы под CC BY-SA, что означает, что вы можете свободно распространять, модифицировать, создавать производные произведения с условием, что вы укажите меня автором оригинального произведения и, в том случае, если вы хотите распространять производную работу, она должна быть под той же лицензией.
Метод молекулярной динамики - метод, в котором временная эволюция системы взаимодействующих атомов или частиц отслеживается интегрированием их уравнений движения (Wikipedia). На практике это означает примерно следующее
Давайте рассмотрим эти положения подробнее.
Надеюсь, вы знаете что такое ньютоновская механика, это что то, строящееся на втором законе Ньютона:
\[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\); воспользовавшись аппаратом статистической физики, мы можем описать распределение частиц по скоростям (например, Максвелловское распределение). Таким образом, можно случайным образом задать частице такую скорость, чтобы температура всей системы была равна заданной. Далее, рассчитать ускорение частицы в нулевой момент времени и, в дальнейшем, собственно координаты частиц в момент времени \(t - \Delta t\).
inb4: Следует понимать, что указанный алгоритм имеет свои недостатки. Например, поскольку скорости явно не учитываются при интегрировании уравнений движения, их расчет имеет некоторую погрешность, что приводит к значительно амплитуде изменения скорости одной частицы со временем (и, как следствие, погрешность при определении кинетической энергии). Также, указанный алгоритм довольно упрощен, для уточнения конкретных реализаций следует читать исходный код (если доступен) соответствующих программных пакетов или научные статьи.
Внимательный (лол) читатель мог обратить внимание, что, даже если рассматривать относительно простые случаи потенциальных энергий - энергии гармонических осцилляторов - они содержат некоторые магические параметры - \(k\) и \(r_0\). Если их физический смысл - я надеюсь - не вызывает вопросов - \(k\) - это что-то похожее на жесткость пружины, а \(r_0\) - расстояние в точке минимума - то вопрос, откуда их брать, может вызвать некоторые сложности.
Для их определения, обычно решается обратная задача - исходя из известной энергии рассчитываются вторые производные (они же \(k\)), а \(r_0\) рассчитывается исходя из минимума потенциальной энергии. В качестве опорных значений для расчетов энергии используются неэмерические методы (в противовес эмпирической классической молекулярной механике), поголовно состоящие из квантово-химических методов.
Совокупность же параметров, рассчитанных в рамках одного метода с усреднением параметров по определенному принципу (равно как и подгон значений под определенные вполне материальные физические величины) в дальнейшем мы будем называть силовым полем. Вопрос выбора того или иного силового поля, а также принципы, которые лежат в основе их определения, мы рассмотрим позже.
В методах статистической физики используются вполне определенные ансабмли. У них есть хитрые названия, которые я использовать не буду, потому что они не очень понятны. Принцип их образования примерно следующий: берутся определенные физические величины, которые нам известны в заданной конфигурации, а остальные физические величины рассчитываются исходя из определенных.
В рамках данных статей я буду рассматривать только те системы, в которых число частиц постоянно. Некоторые величины, например, энергия, довольно трудно зафиксировать, поэтому они также не используются при моделировании. По факту, как правило, фиксируется две из трех величин: давление, температура и объем:
В качестве первого приближения можно считать, что объем можно хитрым образом выразить через давление и температуру в системе (уравнение Менделеева-Клапейрона). Тогда задача сводится к - так или иначе - заданию определенного давления и температуры. Для поддержания определенной температуры в системе используется алгоритм, который описывает термостат (обмен температуры системы с окружающей средой). Для давления, соответственно, баростат - “обмен” давления с окружающей средой.
Под траекторией мы будем подразумевать зависимость координат частиц от времени. Также, в виду лени, сюда мы будем относить и все производные физические величины, которые мы можем рассчитать напрямую в ходе моделирования - скорости, ускорения, температура, давление, различные энергии и т.д.