7 января 2011 в 19:59

Программирование молекулярной динамики из песочницы

Доброе время суток хабражители и хабражительницы. Сегодня хотел бы поделиться с вами своими попытками в программировании физических процессов. А конкретнее – попытке углубиться до молекулярных масштабов. Тема для разговора под хабракатом – молекулярная динамика.

Что есть МД?


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

Математика


Тут все очень просто. Взаимодействия между частицами (в моем конкретном случае атомы) будут считаться с помощью классической физики:

Вот эта, на первый взгляд, запутанная формула – на деле полный аналог школьного F=ma.
(на деле «запутанная» формула длиннее, и выглядит менее эстетично – диссипативную составляющую я намеренно забыл)

Математически, моделирование методом частиц, представляет из себя решение задачи Коши для уравнений выше.
Вычисление, всего этого дела, я поручил алгоритму Верле.

Имхо, он оптимален по точности и скорости. Одно но – ему надо знать два предыдущих положения частицы! Поэтому первый шаг оставим неточному Эйлеру.

Без потенциала никуда


Ф( r ) в формуле выше, определяется через потенциал. Вообще, потенциал тут это царь и бог – именно его влияние самое значительное. Я выбрал довольно популярный — знакомьтесь, Леннард-Джонс:


А вот так он выглядит, здесь а это межатомное расстояние, f — сила взаимодействия, r — расстояние между частицами.

Хватит формул


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

Давайте посмотрим на откольный удар (1083 частицы, ударник «падает» на мишень)











Вот именно этот момент, пожалуй, самый приятный – играться со свежей моделью. Удручала лишь низкая скорость расчета — пришло время оптимизаций.

Ускоряемся


В первую очередь, пришлось вернуться к физике. Лишь один взгляд на потенциал давал понять – при радиусе больше 2a – можно и не считать. Т.е. проще всего тупо для каждой частицы считать взаимодействие лишь с частицами попадающими в сферу радиусом а2 с центром в частице для которой считаем. Одно условие – «если» и расчет ускоряется в разы (зависит от кол-ва частиц и области).

Наблюдать за полной загрузкой одного ядра на моем стационарнике – унылое занятие. Пора бы распараллелить процесс. Тем более именно мд отлично параллелится. Казалось бы все просто – задай соответствие: частица – процессор и считай. Практика показала – лучше паралелить, разбив пространство на сетку и вот почему.
Создав сетку с размером ячейки равным а, мы убиваем сразу двух зайцев – решаем проблему с областью расчета, где учитывается взаимодействие и делаем более прозрачным и простым обмен между процессорами. Вот как это выглядит на примере двух процессоров:


Здесь стоит обмолвится, что любая модель «гонится» еще и различными «хаками» в зависимости от начальных условий. Например, если нам надо смоделировать нечто бесконечное (или огромное) и при этом одинаковой структуры в одном или нескольких направлениях, можно применять периодические граничные условия. Я применял это для моделирования трения двух бесконечных поверхностей.
Смысл тут, в том, что вокруг расчетной области строятся ее «образы», с актуальным положением частиц. И частицы «реальной» области взаимодействуют с частицами в «образе». А если частица пересекает границу расчетной области, она появляется с другой стороны (ака туннель в пакмане). Ниже представлен пример периодических условий для x и y.



А вот для моделирования трения, я использовал периодику только по одной оси:








Вывод



Получилась довольно интересная и наглядная модель, которой вполне можно как играться так и попробовать использовать ее для реального расчета (достоверность результатов очень приличная). По хорошему приложение стоило бы заточить под Cuda или OpenCL — скорость расчета МД на видеокартах, при грамотном распараллеливании, просто запредельная! Впрочем делая ставку на свою лень — скоро это врядли случится.

P.S.: По большей части опирался я лишь на одну книжку: Метод частиц и его использование в механике деформируемого твердого дела. Кривцов А.М, Кривцова Н.В
+71
6716
65
Dare 11,5

Комментарии (32)

+1
GetWindowsDirectory #
А поиграться дадите?
+2
Dare #
В принципе без проблем, если есть желание покопаться в исходниках.

Дело в том, что пока я задаю положения частиц тупо в программном коде. В некоторых случаях это удобно — например «темную» поверхность я задавал используя модуль косинуса (чтобы сделать зубчики).
+3
Bkmz #
Сделайте комментарии для «тупых», «Скорость здесь», «масса здесь».
Будет здорово.)
+2
Dare #
Визуальный редактор уже в альфе.

Там можно будет задавать типы решеток и все такое.
Просто, если тупо натыкать частиц — с вероятностью 99.9% они разлетятся нафиг или будете просто наблюдать газ, что совсем не интересно (надо рисовать устойчивые решетки для твердых тел).
+14
miguello #
Есть куча доступных кодов как для MD, так и для много чего другого типа DFT и так далее. Есть MD, специализированные под разные вещи. Я например, в свое время слегко модифицировал под свои цели (смотрел, как ведет себя углеродная нанотрубка под давлением) код Бреннера: www.eng.fsu.edu/~dommelen/research/nano/brenner/

Большая часть различных кодов написана на Фортране, но есть и на C.
вот пример того, что получилось у меня:




Это, соответственно, компрессия, растяжение, и сгинание (извините, если не совсем по-русски, короче, compression, tension and bending)…

Делалось это кстати все в рамках курса Computational Physics.
+1
miguello #
Упс, сорри за первую картинку, она не должна быть такой большой.
0
ZimM #
Ваш коммент почти на 10 Мб завесил )
0
miguello #
Сорри, картинки не я делал, я делал видео… Ниже отписался… :( Еще раз сорри :)
0
spmbt #
Оптимизировал для Ваc картинки, получилось примерно вдвоее легче. Подмените, пожалуйста, на хостинге исходные файлы теми, что лежат по адресу — читать народу будет легче.

ifolder.ru/21199625
0
miguello #
Сорри, сорри, сорри, картинки не у меня лежат, а у преподавательницы по курсу, я хотлинкнул, я ей сдавал все в виде видео, я так понял, они гифки сделала. Насчет визуального размера — у нее смотрелся нормально, насчет веса — я даже и не ожидал, что можно 10 мегов таким образом набрать. Удалить комментарий тоже не могу, так что сорри, сорри, сорри…
+1
Ember #
А как насчет задачи смоделировать формирование солнца и планет из хаотичной звездной пыли? Берем N-частиц, задаем гравитационное взаимодействие, запускаем и смотрим смогут ли они из хаотичного газа слипнуться в сгусток солнца и планеты вокруг и как будут двигаться?
+1
gribozavr #
Памяти не хватит. Даже если для каждого атома потребуется 1 байт, то чтобы представить 1 моль вещества потребуется 6*10^23/10^9 = 6*10^14 Гб памяти.
+1
Milfgard #
Рассчитывать большими кластерами как метачастицами, фрактализовать?
+2
gribozavr #
Это уже совершенно другой уровень сложности как физических и математических выкладок, так и программирования.
+1
Dare #
Я только о моделировании макроразмеров мечтаю, а вы уже солнечные системы!

Интересная задача для такого метода, но представте какое же количество частиц газа необходимо будет для моделирования целой системы!

Если же мы возьмем абстрактый «булыжник» размером в несколько км как частицу(чтобы сделать реальным время расчета) — тут потребуется уже совсем другой потенциал. Его еще надо придумать.
+2
SeLarin #
Для моделирования макроразмеров надо уже к конечным элементам обращаться… %)
0
Bas1l #
Это лучше стандартными методами гидродинамики обсчитывать. Уравнения Навье-Стокса, метод конечных элементов, и тп. Другой уровень абстракции, так сказать.
0
Nashev #
Если всё это на молекулярном уровне обсчитать, надо бы будет потом приглядеться — не получилось ли на поверхностях некоторых планет тонкой плёнки молекул, скомбинировавшихся в жизнь и цивилизации… %)
+4
SeLarin #
В простейшем варианте MD, действительно, хорошо ускоряется на GPU, но только до тех пор пока на арену не выходят дальнодействующие взаимодействия, например, электростатические. Корректный учет зарядов в системах приближенных к реальным — это огромная головная боль и существующие алгоритмы, позволяющие делать это не так уж и хорошо параллелятся на той же CUDA.

А вообще в плане практического моделирования такого рода систем я бы посоветовал посмотреть на пакет ESPResSo. Он многофункциональный, но если пошаманить опциями при сборке и отключить все ненужное, то можно получить довольно шустрое создание, параллельное из коробки (через MPI). У него есть только один минус: все параметры, сам цикл моделирования и обработки надо будет писать на Tcl…
+2
youROCK #
Автор, потенциал равен интегральчику от f®, ну или сила равна производной от f®, а не тупо 1/r * f® ;). Ну и да, ИМХО, можно было бы использовать существующие решения для МД, например LAMMPS, а писать самому только визуализацию :). Потому что с интегрированием уравнений движения не всё так просто, надо использовать правильные схемы, которые являются хотя бы консервативными, а не просто Эйлера (или кого-нибудь другого :)). А в остальном всё замечательно :)
0
youROCK #
Хабрапарсеру привет… f ( r ), а не f®:
0
youROCK #
А вообще, формула немного странная у Вас используется, ИМХО

Всё же вроде как очень просто —


И не надо никаких дополнительных и сбивающих с толку формул :)
0
SeLarin #
А какая формула странная? Если потенциалы, то там, как я понимаю, написаны формулы для конкретного случая, взятые из книжки (надо смотреть книжку, чтобы точнее сказать). Если «странная» формула, описывающая алгоритм Верле, то это вроде как классика…
0
youROCK #
Самая первая, которая про второй закон Ньютона :). ИМХО, она написана не совсем корректно, и я привел тот вариант, который, ИМХО, является более корректным.
0
Dare #
Если честно — я просто заскриншотил кривцова, по которому писал.
Да, там чуток заморочено
0
CGS #
… диссипативную составляющую ...

Давно меня так не торкало…
+1
mr_idiot #
У меня диплом по молекулярной динамике был. В моей проге можно создавать различные решетки, задавать потенциалы и делать другие штуки. Могу выложить.
0
mr_idiot #
ifolder.ru/21204653 — вот такая вот программа. На ней я тоже всякие исследования делал. Тоже статью думал написать, но было лень.
0
Dare #
Очень интересно! Может быть все же напишете статью?)
Считает программа у вас довольно шустро, хотя судя по всему не параллелена.
+1
Bas1l #
А вы не задумывались насчет GROMACS?

Кстати, по молекулярной динамике есть хорошая книжка Computer simulation of liquids (и <a href=«www.filestube.com/abb88c84ad44045503ea,g/Computer-Simulation-of-Liquids.html»бесплатно). Она 1989 года, но нам ее рекомендовал в универе бывший препод, который сейчас профессионально занимается MD в Германии. Осторожно: ОЧЕНЬ много математики (даже кватернионы есть:-)).

Кстати, молекулярная динамика--это не только расчет траектории движения одной системы в фазовом пространстве, но и расчет распределения систем по состояниям фазового пространства при заданных внешних условиях. И там совсем другие методы.
0
Dare #
Нет, о пакетах не задумывался вообще. Хотелось для начала, самому написать.

Книжку уже качаю — после сессии взгляну обязательно, сейчас мозг на другое переключен.

Спасибо!
–1
auraz #
интересная статья, спасибо

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.