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

Что есть МД?

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

Математика

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

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

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

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

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

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


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

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

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

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











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

Ускоряемся

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

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


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



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

Вывод

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

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