Событийный подход в молекулярной динамике

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

Сегодня метод молекулярной динамики является одним из доминирующих методов компьютерного моделирования физических процессов в мировой компьютерной науке [1]. Дополненный методами теоретической физики, теории устойчивости, случайных процессов и др., метод молекулярной динамики является основным инструментом исследований систем с большим числом степеней свободы. В настоящее время ведется активное развитие методов квантовой молекулярной динамики.

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

Весьма эффективным является событийный подход в молекулярной динамике. Его применяют для много компонентных систем, составляющие которых развиваются независимо. В [2] рассматривается система, состоящая из N твердых сфер в прямоугольной области с периодическими граничными условиями. В этом алгоритме рассматриваются два основных объекта данных: состояния частиц и список событий. Все события записываются в очередь по приоритету. Сложность извлечения следующего события из такой очереди - , вставки нового события - . Пространство моделирования делят на ячейки одинакового размера.

1. СОБЫТИЙНЫЙ ПОДХОД В МОЛЕКУЛЯРНОЙ ДИНАМИКЕ

В событийной молекулярной динамике асинхронно обрабатывается последовательность событий. Прямой, но упрощенный подход обрабатывает все частицы при каждом событии. В более сложном алгоритме обрабатываются только те частицы, которые задействованы в событии. Этот подход был успешно реализован для решения многих задач, среди которых гранулированное топливо [7], привязные мембраны [7], модели полей сражений.

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

Единственная причина, по которой событийная молекулярная динамика не стала также популярна, как традиционная молекулярная динамика, управляемая временем, это то, что распараллеливание является крайне сложной процедурой. Парадоксальность этого задания заключается в алгоритмическом распараллеливании физически непараллельной динамики. Тем не менее, некоторые идеи и основные положения были предложены в работе [8]. Они особенно подходят для компьютеров с разделяемой памятью, но могут также быть использованы для архитектур с разделенной памятью.

Рис. 1 - Схема основной части алгоритма

,

где - время состояния, - текущее время, - время события.

Рис. 2 - Структура данных для N частиц

Состояние частицы описано ее физическими свойствами, такими как координаты, скорость и др., а после обработки события времени этого события - типа события, номера ячейки, глобального номера частицы, который необходим для идентификации частиц при обработке разных событий.

Список событий ставит в соответствие каждой частице событие в будущем. Единица данных списка событий состоит из времени события, типа события, соучастника события, и копии счетчика событий для соучастника.

Все события записываются в очередь по приоритету, которая обычно реализуется как имплицитная куча [8, 9], но применение других структур данных с подобными характеристиками также возможно [10]. Куча представляет собой двоичное дерево со следующим правилом упорядочивания: каждый родительский узел имеет более высокий приоритет, чем его потомки; потомки между собой не упорядочены. Таким образом, корень всегда имеет высший приоритет, то есть самое раннее время события. Получение следующего события из очереди с приоритетом имеет вычислительную сложность . Вставка нового события в очередь по приоритету имеет сложность .

Такая организация данных является более эффективной по сравнению с алгоритмами в [7], так как здесь нет промежуточного хранения данных, хотя основные идеи [7] все еще имеют силу.

,

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

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

Более ранний алгоритм проверяет при каждом столкновении, становится ли соучастник события недействительным и, если да, помечает этого партнера. Эта стратегия является менее эффективной, так как иногда одна и та же частица может быть определена как недействительная несколько раз. Кроме того, в параллельном алгоритме, возможно, будет необходимым дополнительное сообщение между различными процессами.

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

Две частицы 1 и 2 столкнуться в момент времени :

,

где , - это относительные скорости и положения частиц в момент времени , а - радиусы частиц. Если мнимо или меньше , частицы не столкнутся.

Если алгоритм должен осуществлять проверку столкновений со всеми частицами, то его скорость выполнения будет очень низкой для большого числа частиц. Поэтому мы делим пространство модели на С ячеек равного размера. Каждая частица принадлежит ячейке, в которой лежит ее центр. Если размер ячейки больше, чем максимальное расстояние между двумя частицами при взаимодействии, то есть в случае твердых сфер превышает диаметр частиц, вычисление событий должно происходить только тогда, когда частицы лежат в одной или соседних ячейках (одна квадратная ячейка имеет 8 соседних в двумерном случае и 26 в трехмерном).

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

В среднем для каждой ячейки в соседних для нее находится частиц. Если , это количество сильно возрастает и после каждого столкновения приходится пересчитывать большое количество возможных событий. С другой стороны, если , то каждой частице приходится пересекать большое количество граничных ячеек в промежутке между столкновениями двух частиц, что приводит к большему количеству обрабатываемых событий для полного моделирования. Эти два замечания противопоставлены друг другу. Первое пропорционально , а второе - . На рис. 3 (слева) показано, что существует нечеткий минимум для времени симуляции при низкой плотности при оптимальном количестве ячеек в двумерном случае и в трехмерном. Программа может выбрать оптимальное количество ячеек до начала вычислений. Заметим, что в пределе для высоких плотностей, особенно в трехмерном случае, оптимальное количество ячеек не может быть достигнуто, так как размеры ячейки должны быть больше, чем диаметры частиц (рис. 3 справа). В этом случае должно быть выбрано как можно большим.

Рис. 3 - График зависимости времени симуляции от количества ячеек

В примере (рис. 3) использованы частиц (двумерный случай) и частиц (трехмерный случай) с относительными объемами (двумерный), (трехмерный случай на левом рисунке) и (трехмерный случай на правом рисунке).

2. МОДЕЛЬ ТЕЧЕНИЯ ГАЗА В СРЕДЕ С ФИЛЬТРАМИ

Рассмотрим двумерный случай течения жидкости в некоторой прямоугольной области. Определим, что вектор движения жидкости направлен слева направо. Будем считать, что рассматриваемая область является замкнутой, то есть молекулы жидкости, вылетая из области на правом срезе, передаются на вход (левый срез). Появление на входе частиц осуществляется по принципу очереди «Первый зашел, первый вышел». Границы области будем считать отражающими. В рабочей области расположены фильтры. Каждый фильтр - это набор одинаковых прямоугольников с отражающими границами, расположенных параллельно горизонтальным границам области. Расстояния между прямоугольниками постоянны. Пример области с фильтрами показан на рис. 4.

Рис. 4 - Рабочая область с фильтрами

Кроме основного вещества в области также движутся включения - непроницаемые тела формы, близкой к прямоугольной, обладающими свойствами упругости и гибкости в той или иной степени. Эти включения перемещаются слева направо вдоль рабочей области, так как на них действуют силы со стороны основного потока. Они могут быть остановлены при столкновении с фильтром. При большом количестве включений движение основного потока может прекратиться в силу закупорки фильтров.

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

Частицы жидкости при столкновении упруго отражаются. В зависимости от типа столкновения будем различать внутренние и внешние столкновения.

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

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

На этапе инициализации множество частиц, границ ячеек и области генерируется автоматически.

Согласно рассмотренному в первом разделе алгоритму (рис. 1.1) на каждом шаге мы получаем из хипа («кучи») время ближайшего события и номера частиц, для которых оно должно произойти. Хип - это дерево, которое формируется из локальных времен всех частиц по следующему правилу: значение отцовского узла не превышает минимального времени до ближайшего события рассматриваемых частиц. Его инициализация происходит до первого прохождения алгоритма и в дальнейшем алгоритм лишь добавляет новые значения, убирая старые. В случае если количество модельных частиц отлично от некоторой степени двойки, то удобно дополнить дерево фиктивными значениями, заведомо большими любого из времен наступления события. В таком случае они не повлияют на результат.

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

Работа алгоритма продолжается до тех пор, пока не будет остановлена.

Выводы

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

В данной работе рассмотрен процесс моделирования потока жидкости в замкнутой области с фильтрами. Для большей наглядности область была представлена прямоугольной формы. Для моделирования был использован событийный подход, основанный на пересчете столкновений (событий) для модельных частиц и обработке лишь тех, для которых момент наступления событий - минимален.

Список использованной литературы

1. Allen P., Introduction to Molecular Dynamics Simulation - Computational Soft Matter: From Synthetic Polymers To Proteins, Lecture Notes, 2004 - p. 1-24.

2. Miller S., Luding S, Event-driven molecular dynamics in parallel / S. Miller, S. Luding - Journal of Computational Physics, 193, 2004 - p. 306-316.

3. Kuksin A., Morozov I., Norman G., Stegailov V., Valuev I., Standard of molecular dynamics modeling and simulation of relaxation - Molecular simulation, 2005. V. 31. - p. 1005-1017.

4. De Michele C., Scala A., Schilling R., Sciortino F., Molecular correlation functions for uniaxial ellipsoids in the isotropic state, - J. Phys. Chem., 124:104509, 2004.

5. Scala A, De Michele C., Voightmann Th., Event-driven Brownian dynamics for hard spheres - ArXiv Condensed Matter e-prints, 2006.

6. S. Luding, Computer Physics Communications 147, 134, 2002.

7. B. D. Lubachevsky, Int. J. Comput. Simul. 2, 372, 1992.

8. D. E. Knuth, The Art of Computer Programming - Addison-Wesley, Reading, MA, 1968.

9. Чернышев Ю.К. Событийное программирование. Применение к решению некоторых задач физики / Ю.К. Чернышев // НАУ им. Н.Е. Жуковского «ХАИ». -- Х., 2008.