azTotMD

Методы численного моделирования молекулярных систем, такие как молекулярная динамика, рассматривают эти системы как механические (что-то вроде набора шариков на пружинках). Однако, в отличие от механических систем, для молекулярных существует понятие температура. Вещество не может существовать без температуры, а температура – без вещества (на счет последней части утверждения есть и другое мнение). Из опыта мы знаем, что очень многие свойства вещества кардинально зависят от температуры, и, естественно, что её надо как-то учитывать. Для поддержания температуры в молекулярной динамике используются специальные алгоритмы-«термостаты». Наиболее известные среди них это термостаты Андерсена, Берендсена и Нозе-Гувера. Все они основаны на молекулярно-кинетической теории газов, где температура есть просто величина пропорциональная среднекинетической энергии молекул. Соответственно, работа данных термостатов осуществляется путём умножения скоростей частиц на некоторую величину.

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

Например, вода:



очевидно, что водород с кислородом внутри одной молекулы взаимодействуют совсем по-другому, нежели тот же кислород с водородом соседней молекулы. Таким образом, разделяют ВНУТРИмолекулярное (intramolecular) и МЕЖмолекулярное (intermolecular) взаимодействие. Межмолекулярное взаимодействие можно задать короткодействующими и Кулоновскими парными потенциалами, о которых речь шла в предыдущих статьях. Здесь же сконцентрируемся на внутримолекулярном.

В предыдущей статье мы обсудили основу метода молекулярной динамики, в том числе вычисление энергии и сил взаимодействия между частицами с заданными парными потенциалами. А что, если частицы обладают некоторым электрическим зарядом? Например, в том случае, если мы моделируем кристалл поваренной соли, состоящий из ионов Na⁺ и Cl^-. Или водный раствор, содержащий те или иные ионы. В этом случае, кроме парных потенциалов типа Леннарда-Джонса между ионами действуют силы электростатического взаимодействия, т.е. закон Кулона. Энергия такого взаимодействия для пары частиц i-j равна:

$$

где q – заряд частицы, r_ij – расстояние между частицами, С – некоторая постоянная, зависящая от выбора единиц измерения. В системе СИ это — $$, в СГС — 1, в моей программе (где энергия выражена в электронвольтах, расстояние в ангстремах, а заряд в элементарных зарядах) C примерно равно 14.3996.



Ну и что, скажете вы? Просто добавим соответствующее слагаемое в парный потенциал и готово. Однако, чаще всего в МД моделировании используют периодические граничные условия, т.е. моделируемая система со всех сторон окружена бесконечным количеством её виртуальных копий. В этом случае каждый виртуальный образ нашей системы будет взаимодействовать со всеми заряженными частицами внутри системы по закону Кулона. А поскольку Кулоновское взаимодействие убывает с расстоянием очень слабо (как 1/r), то отмахнуться от него так просто нельзя, сказав, что с такого-то расстояния мы его не вычисляем. Ряд вида 1/x расходится, т.е. его сумма, в принципе, может расти до бесконечности. И что же теперь, миску супа не солить? Убьёт электричеством?

Цель данной статьи – поднять вопросы распараллеливания кода программы для численного моделирования методом молекулярной динамики (МД) с помощью технологии CUDA. Зачем это вообще нужно, ведь уже существуют программные пакеты по МД, работающие в том числе и на CUDA? Дело в том, что я развиваю свою собственную концепцию «непостоянного поля сил» (non-constant force field), которая не реализована в существующих МД-программах.

Переделывать чужой код под эти нужды – довольно неблагодарное занятие, поэтому я взялся перенести уже написанный свой последовательный код и заодно поделится некоторыми размышлениями. Кроме того, это ответ на часто мелькающий здесь комментарий к статьям по CUDA, вроде этого .

Итак, что же такое молекулярная динамика? На Хабре уже есть несколько постов на эту тему, например здесь или вот здесь. Кратко, МД – это метод, позволяющий моделировать движение множества частиц (в том числе атомов, ионов, молекул) и рассчитывать коллективные свойства системы, зависящие от этого движения. Как это работает? Допустим для множества из N частиц заданы некоторые начальные координаты, скорости, массы и (главное!) законы взаимодействия между ними. Изменяем координаты согласно скоростям. На основе законов взаимодействия вычисляем силы, действующие между частицами. Раз знаем силу и массу – знаем ускорение. Поправляем скорость с учетом ускорения. И снова переходим к изменению координат. И так повторяем тысячи раз, пока не надоест не наберем достаточную статистику.