Содержание к диссертации
Введение
Глава 1. Моделирование задач молекулярной динамики 15
1.1 Алгоритмы молекулярной динамики 16
1.2 Граничные условия 17
1.3 Начальные условия 20
Начальные координаты
Начальные скорости
1.4 Потенциалы взаимодействия 28
1.5 Термостаты 36
1.6 Численные методы 37
Глава 2. Постановка задачи установления термодинамического равновесия 44
2.1 Постановка задачи для аргона 45
2.1.1 Граничные условия 46
2.1.2 Начальные условия 47
2.1.3 Потенциалы 49
2.1.4 Численная схема решения 52
2.2 Постановка задачи для алюминия 54
2.2.1 Граничные условия 56
2.2.2 Начальные условия 56
2.2.3 Потенциалы 60
2.2.4 Алгоритм нагрева (термостат) 63
2.2.5 Численная схема решения 64
Глава 3. Программная реализация алгоритма и параллельные вычисления 67
3.1 Тонкости программной реализации 68
3.2 Параллельные вычисления 72
Глава 4. Физическая задача, анализ данных, результаты вычислений 84
4.1 Задача для аргона 84
4.2 Задача для алюминия 89
Заключение 98
Список литературы 100
- Потенциалы взаимодействия
- Численная схема решения
- Алгоритм нагрева (термостат)
- Тонкости программной реализации
Введение к работе
Актуальность. Современная вычислительная техника дает возможность с помощью компьютерного моделирования решать недоступные для аналитического решения задачи на микроскопическом уровне. Компьютерная молекулярная динамика является одним из наиболее мощных вычислительных методов, эффективно применяемых для моделирования физических, химических и биологических систем. В молекулярной динамике (МД) исследуемая система представляется совокупностью взаимодействующих атомов или молекул, движущихся по законам классической механики. Благодаря развитию и применению компьютерных технологий и графических методов анализа эффективность применяемых методов МД-моделирования неуклонно возрастает. Моделирование реальных физических систем, например, кристаллов или огромных биологических молекул на базе методов МД представляет собой определяющее перспективное направление в ближайшем будущем. Необходимым условием для расчета задач моделирования динамических процессов, решаемых с помощью метода молекулярной динамики, является равновесное макроскопическое состояние системы при выбранных тепловых условиях в начальный момент времени. Получение характеристик состояния термодинамического равновесия системы является отдельной задачей, и чаще всего их получают экспериментально, что является дорогостоящим, трудоемким и длительным процессом.
Одним из основных этапов решения задач в молекулярной динамике является выбор модели взаимодействия атомов системы, которая позволит получать в поставленных условиях реалистичные данные, наиболее схожие с результатами эксперимента. Взаимодействие частиц описывается с помощью потенциала, который определяется исходя из свойств моделируемого вещества, его агрегатного и термического состояний. Несмотря на существование большого количества различных видов потенциалов, ограничения в подборе параметров и их неспособности воспроизвести полный набор физических характеристик порождают существенные сложности в выборе подходящего потенциала при моделировании процессов в узконаправленных задачах.
На каждом шаге моделирования решается система обыкновенных дифференциальных уравнений второго порядка, соответствующая второму закону Ньютона и описывающая движение частиц молекулярно-динамической системы. При решении задачи необходимо подобрать оптимальную численную схему, аппроксимирующую уравнения движения. Выбор схемы численного интегрирования обуславливается требованиями устойчивости и минимальности затрат машинного времени при расчете точек фазовой траектории.
Важным моментом является выбор временного шага интегрирования уравнений движения. Следует учитывать, что достаточно большая его величина приводит к несохранению полной энергии и неустойчивым численным решениям, то есть к таким, которые все больше отклоняются с течением времени от истинного решения. С другой стороны, размер шага интегрирования определяет реальное время моделирования системы. Существенно малый временной шаг может привести к тому, что за время компьютерного эксперимента система не достигнет равновесия.
Одними из основных инструментов моделирования задач МД в России и мире являются пакеты прикладных программ, ориентированные на высокопроизводительные вычислительные системы. Ввиду развития суперкомпьютеров, построенных на основе графических плат, встает вопрос о переносе пакетов прикладных программ на новые платформы для повышения эффективности их использования и расширения диапазона решаемых задач за счет роста производительности вычислительных мощностей. По этой причине актуальным является развитие собственных компьютерных разработок -разработок программных комплексов, позволяющих моделировать задачи молекулярной динамики и являющихся более доступными в использовании.
Следует отметить, что в настоящее время все больший интерес у специалистов по математическому моделированию вызывают гетерогенные вычислительные системы, основанные на графических ускорителях; алгоритмы молекулярной динамики являются одними из алгоритмов, эффективно реализующихся на подобных архитектурах.
Вышесказанное приводит к необходимости разработки новых математических моделей, новых алгоритмов, параллельных программ и методов анализа результатов для исследования задач методом молекулярной динамики.
Цели работы.
Моделирование процесса установления термодинамического равновесия систем веществ, имеющих разные структуры, методами молекулярной динамики.
Исследование зависимости структуры и характеристик рассматриваемого вещества от геометрии расчетного объема, количества частиц и температуры.
Разработка алгоритмов моделирования и анализа результатов молекулярных расчетов для систем атомов в разных агрегатных состояниях.
Создание программной реализации с использованием высокопроизводительных технологий.
Научная новизна. Предложены и реализованы алгоритмы по моделированию процесса установления термодинамического равновесия методами молекулярной динамики систем веществ, имеющих разные структуры (газы, металлы).
Теоретическая и практическая значимость.
Создана оригинальная молекулярно-динамическая программа, позволяющая изучать зависимость структуры и характеристик исследуемого вещества от геометрии расчетного объема, разного количества частиц и от температуры.
Получены термодинамические характеристики для систем аргонового газа в объемной геометрии и алюминия в геометрии пластины.
Предложен и реализован адаптивный алгоритм подбора шага интегрирования, позволяющий обеспечивать контроль над изменением полной энергии системы в течение всего времени расчета.
Для определения значения средних температур в данной работе конструируются методы анализа, основанные на изучении экспоненциального показателя равновесного распределения.
Проведен анализ временного параметра алгоритма теплового контроля (термостата) для эффективного получения желаемой температуры в образце.
Разработаны и реализованы алгоритмы параллельных расчетов на гибридных вычислительных комплексах при моделировании установления термодинамического равновесия систем веществ в разных агрегатных состояниях. Параллельная реализация адаптирована для применения технологии CUDA с учетом эффективного распределения данных.
Апробация работы. Основные результаты диссертационной работы докладывались и обсуждались на следующих конференции и семинарах:
XVIII-й Международной научной конференции студентов, аспирантов и молодых ученых «Ломоносов - 2011» в секции «Вычислительная математика и кибернетика», МГУ имени М. В. Ломоносова (Москва, апрель 2011);
научном семинаре «Математическое моделирование» под руководством профессора В. Ф. Тишкина и профессора А. А. Кулешова в ИПМ имени М. В. Келдыша РАН (Москва, май 2011);
научном семинаре кафедры вычислительных методов факультета ВМК МГУ имени М. В. Ломоносова под руководством профессора А. В. Гулина (Москва, май 2011);
научном семинаре кафедры информатики МФТИ под руководством профессора И. Б. Петрова (Москва, сентябрь 2011).
Публикации. Основные результаты диссертации опубликованы в трех работах, две из которых в изданиях, рекомендованных ВАК [1, 3].
Структура и объем диссертации. Диссертация состоит из оглавления, введения, четырех глав, заключения и списка литературы. Полный объем диссертации составляет 114 страниц, включая 21 иллюстрацию и 2 таблицы, пронумерованные по главам. Список литературы содержит 112 наименований.
Потенциалы взаимодействия
Потенциальная энергия - функция, зависящая от координат рассматриваемых частиц и описывающая взаимодействие между частицами системы. На данный момент времени известно множество форм потенциалов межатомного взаимодействия [11, 30, 53, 87, ПО]. По учету действующих величин они разделяются на потенциалы парного взаимодействия и потенциалы многоатомного взаимодействия.
Парные потенциалы рассматривают взаимодействие двух частиц, сама функция зависит только от координат рассматриваемых атомов. Основным свойством парных потенциалов является отталкивание при сближении и притягивание частиц при удалении, причем при значительном удалении значение потенциальной энергии становится пренебрежимо мало, то есть потенциальная функция "стремится" к нулю.
Для многоатомных потенциалов [18, 69] характерно рассмотрение более двух частиц во взаимодействии, учитывание взаимной ориентации молекул, а именно углов между связями. Многочастичные потенциалы [16] используются в моделях большинства неплотноупакованных решеток и многоатомных молекул, поскольку расстояния между дальнейшими соседними атомами оказываются на неустойчивом участке парного потенциала взаимодействия и потому парные потенциалы не способны корректно описать взаимосвязь в таких системах. Многоатомные потенциалы получили большое распространение при описании молекулярных систем, но такая модель часто оказывается сугубо эмпирической и требующей подбора большого числа констант, справедливых только для данной определенной задачи с конкретным соединением. Главным недостатком многочастичных потенциалов является то, что при диссоциации молекул и разрушении кристаллических решеток они теряют физический смысл, а-следовательно сфера их возможного применения весьма ограничена.
Выбор определенного вида потенциала взаимодействия основывается на сравнении механических свойств компьютерной модели потенциала и реального материала.
Модельные потенциалы зависят от ряда своих параметров, которые подбираются таким образом, чтоб вычисленные с его помощью значения физических величин совпадали или были достаточно близки к значениям, определяемым экспериментальным путем. Потенциальная энергия /-ой частицы Ui при использовании в моделировании парного потенциала и(гу) представляется в виде суммы парных взаимодействий: Потенциал Леннарда-Джонса (1.10) является двухпараметрическим парным потенциалом, что ограничивает возможности удовлетворения сочетаниям макроскопических параметров, а также не способствует его использованию при моделировании тел с большими плотностями, где взаимодействие не является парным, так как должно учитываться влияние окружений частиц. Несомненным достоинством его является простота вычислительной формы и то, что он достаточно точно описывает свойства веществ, для которых является возможным парное взаимодействие (например, кристаллических инертных газов). Еще одним часто используемым потенциалом является трехпараметрический парный потенциал Морзе [90]: р — параметр жесткости, от этого параметра зависит поведение потенциала в области гу а. Потенциал Морзе (1.11) является трехпараметрическим, что дает больше возможностей для вариаций макроскопических переменных по сравнению с потенциалом Леннарда-Джонса (1.10). При определенном подборе коэффициентов (а-/? = 6) эти два потенциала становятся близки по характеристикам [25]. Еще одним достоинством потенциала является экспоненциальная форма, которая означает быстрое затухание потенциала на расстоянии, потенциал Леннарда-Джонса имеет степенной вид и поэтому затухает более медленно, но такая форма записи усложняет компьютерный расчет, так как вычисления функции экспоненты и квадратного корня для определения длины вектора расстояния между частицами являются трудоемкими. Парный потенциал Морзе удовлетворительно описывает свойства веществ, имеющих гранецентрированную кубическую решетку [86]. Известно, что для большинства металлов и полупроводников парные потенциалы не могут обеспечить реалистичных значений [47, 67]. В металлах учет только взаимодействия с ближайшим соседом не достаточен из-за коллективизации электронов проводимости. Часто оказывается необходимым учитывать во взаимодействии около 20 ближайших атомов в цепочке. Для корректного описания свойств твердого тела в большинстве случаев используют эмпирические многочастичные потенциалы [51], которые имеют ограничения в виду подборки параметров для узконаправленной задачи и своей неспособности воспроизвести полный набор физических характеристик. По этой причине выбор необходимой потенциальной функции зависит от желаемых исходных данных и от конкретной физической задачи. Примерами таких потенциалов могут быть потенциал погруженного атома [6, 55, 66, 71], потенциал Финниса-Синклера [70] потенциал Терсоффа [103, 104] , потенциал Стиллингера-Вебера [100].
Численная схема решения
Для решения поставленной задачи используем следующий алгоритм: 1) Расчет по формуле (2.18) новых координат частиц г" для момента времени t = t" + At, где п - номер шага. 2) Проверка граничных условий (2.6). 3) Вычисление значения парциальной U"+l (2.11) и суммарной потенциальной энергии системы U" (2.3) частиц и сил F, (2.16), действующих на частицы. 4) Расчет новых скоростей частиц V"+l для момента времени tn+1 по выражению (2.19). 5) Вычисление кинетической энергии системы КЕп+і по формуле (2.4). 6) Вычисление полной энергии системы En+l (2.5) для момента времени tn+i. Этапы 1) - 6) повторяются на каждом шаге по времени. Существует множество способов численного интегрирования уравнений движения [36, 38, 96]. Проблема метода частиц заключается в необходимости интегрирования очень большого числа уравнений. При расчетах основное время уходит на вычисление силы, действующей на данную частицу, что снижает эффективность методов, требующих на каждом шаге многократного вычисления правой части уравнений, и делает привлекательным алгоритм Верле в скоростной форме, для которого значение сил взаимодействия для каждой частицы вычисляются один раз в течение одного шага. С целью слежения за изменением полной энергии системы в работе предлагается следующий алгоритм: характеризующую Вычисляем на каждом шаге величину 5Е относительное изменение полной энергии за время одного шага. Введем ограничение на эту величину сверху и снизу и свяжем ее с изменением временного шага согласно следующему правилу: 3) Если 10 6 5Е 1(Г , то продолжаем счет программы с прежним шагом.
Необходимо отдельно отметить присутствие в знаменателе величины Е". При большом числе частиц суммарная погрешность вычислений разности энергий на новом и предыдущем шагах может принимать большие значения, нормирование погрешности на величину Е" позволяет установить универсальные значения для допустимых границ ее изменения. Это позволяет обеспечивать контроль над изменением энергии в широком диапазоне ее значений.
Алгоритм, составленный по рассматриваемой постановке задачи, лег в основу молекулярно-динамической программы, позволяющей рассчитывать характеристики газов на основе потенциала Леннарда-Джонса. В результате проведенных расчетов была исследована зависимость скорости установления термодинамического равновесия системы, состоящей из атомов аргона, для разного числа частиц.
Установление термодинамического равновесия для системы частиц алюминия при необходимой температуре) При моделировании задачи по установлению термодинамического равновесия твердых материй необходимо выбрать геометрию объекта, то есть определить будет это объемная модель или какая-то другая. Представляет собой как теоретический, так и практический интерес исследование металлов в геометрии пластины. Далее изучаются известные характеристики рассматриваемого материала, по которым выстраивается область, определяется функция взаимодействия частиц в системе, для металлов парные потенциалы не дают корректных данных, используются многоатомные потенциалы, которые учитывают влияние окружение на частицу. Состояние равновесия системы атомов металлов является обязательным условием для расчета задач, связанных с получением свойств металла при критических условиях, таких как изучение процесса абляции под действием фемтосекундного лазера, разрушения структуры металла; характеристик плавления металлов — задач с изменением термических условий. Для задач такого рода, которые являются актуальными в современных исследованиях, образец должен быть подготовлен специальным образом, а именно находиться в состоянии термодинамического равновесия при выбранной температуре, необходимой для начала расчета. Для такой подготовки подбираются в зависимости от условий задачи численные алгоритмы нагрева системы — термостаты. Термостаты являются внешними факторами воздействия и должны учитываться в схеме интегрирования уравнений движения.
Таким образом, математическая постановка задачи состоит в необходимости интегрирования уравнений движения частиц, распределенных в соответствии со структурой материала, с заданным потенциалом взаимодействия при выбранном условии на границах и с учетом воздействия термостата. Рассматривается динамика для алюминиевого образца в геометрии пластины. В такой постановке необходимо выбрать одну из осей координат, по которой модель будет конечной, на две другие оси накладывается условие периодичности, рассмотренное подробно в пункте 2.1.1.
Для моделирования бесконечной пластины вдоль координатных осей х и у наложим условия периодичности в направлении осей с периодами Lx и Ly соответственно, в направлении оси z ограничим систему. Задача решается в области конечного объема расчетной ячейки.
Алгоритм нагрева (термостат)
Математическое моделирование свойств веществ и создание соответствующих комплексов программ являются составной частью приоритетной задачи современной науки.
Потребности новейших технологий привели к исследованию ситуаций, для которых физические эксперименты крайне затруднены или неосуществимы, а чисто теоретический анализ слишком сложен. В настоящее время, благодаря стремительному развитию вычислительной техники, решать такого рода проблемы позволяет математическое моделирование с применением ЭВМ, которое во многих случаях оказывается единственным способом получения детальных количественных сведений о поведении сложных молекулярных систем, известных в природе или, что важнее, еще только планируемых к созданию. Компьютерное моделирование часто выступает в качестве связующего звена между теорией и физическим экспериментом. Это особенно важно тогда, когда в теории не удается отыскать точного аналитического решения сложных систем уравнений, описывающих поведение молекул.
Одним из перспективных и активно развивающихся методов компьютерного моделирования является метод молекулярной динамики, позволяющий определить целый комплекс свойств (структурные, термодинамические, транспортные) и исследовать их взаимосвязи. При этом точность получаемых результатов определяется используемыми математическими моделями и размерностью (числом частиц) моделируемой системы.
В методе молекулярной динамики молекула рассматривается как система взаимодействующих классических частиц. Это означает, что атомы, из которых состоит молекула, считаются шариками, взаимодействующими между собой по законам классической механики. Метод основан на численном интегрировании уравнений движения для системы N частиц.
Время, необходимое для расчета траектории молекул, зависит от размера моделируемой системы и величины шага интегрирования. Чем больше рассматриваемая система, тем точнее будут получаемые результаты, но увеличение размера системы сопровождается значительными тратами компьютерного времени, связанными с возросшим количеством рассматриваемых пар молекул. Чем больше шаг интегрирования, тем меньше шагов потребуется проделать процессору для расчета молекулярной траектории заданной длины, но от величины шага зависит точность проделанных вычислений: чем шаг больше, тем точность меньше.
Уменьшить затраты компьютерного времени, не пренебрегая точностью вычислений, можно, используя периодические граничные условия и оптимизируя вычисления сил взаимодействия, которые представляют собой наиболее затратные части алгоритма расчета траекторий. Применяя периодические граничные условия, реальную бесконечную систему заменяют совокупностью одинаковых конечных систем, каждая из которых расположена в пространственной ячейке, окруженной со всех сторон тождественными ячейками. Молекулы, находящиеся в одной ячейке, взаимодействуют друг с другом и молекулами в соседних ячейках. Взаимодействующими считаются молекулы, расположенные внутри сферы с радиусом гс (радиус сферического отсечения взаимодействий или радиус «обрезания» потенциала). Взаимодействие с молекулами, расположенными вне этой сферы, не учитывается. Внутрь сферы отсечения включаются молекулы, расположенные как в центральной ячейке, так и в соседних ячейках (рис. 3.1). Моделируется динамика лишь одной такой ячейки. Размер ячейки должен быть достаточно большим для исключения возможности краевых эффектов -молекула не должна взаимодействовать сама с собой.
Использование периодических граничных условий при ограничении расстояния взаимодействия радиусом сферического отсечения заключается в том, что для всех пар частиц вычисляются расстояния между ними, а силы воздействия рассчитываются только для пар атомов, расстояние между которыми меньше радиуса «обрезания» потенциала, что значительно снижает количество операций и соответственно время расчета.
Алгоритм нахождения расстояния между двумя частицами тоже является весьма затратным, так как приходится для каждой пары частиц выполнять ряд арифметических операций в виде вычисления квадратного корня из суммы квадратов разностей координат, для того чтобы провести сравнение полученного значения с известной константой гс.
С целью решения данной проблемы и минимизации затрат был реализован метод геометрического разбиения моделируемой области с учетом периодических условий на границах. Область разбивается на «ящики» кубической формы с тремя измерениями, равными радиусу «обрезания». Тогда если рассматриваемая частица находится в определенном «ящике», то считается, что на нее могут действовать атомы из этого же «ящика» либо соседних «ящиков». В двумерном случае приходится проводить расчеты всего для 9 таких участков (рис. 3.2), в трехмерном - для 27 участков, что сокращает компьютерные затраты на порядки, и тем выше порядок сокращения затрат, чем больше размер системы.
Тонкости программной реализации
Представленная компанией NVIDIA программно-аппаратная архитектура для расчетов на графических процессорах CUDA хорошо подходит для решения широкого круга задач с высоким параллелизмом. CUDA работает на большом количестве продуктов NVIDIA, и улучшает модель программирования GPU, значительно упрощая ее и добавляя большое количество возможностей, таких как возможность синхронизации потоков, вычисления с двойной точностью и целочисленные операции.
Возможность создания и использования суперкомпьютеров во всем мире относится к факторам стратегического потенциала научно-технического значения, т.к. прогресс любой развитой страны в современном мире невозможен без компьютеризации во всех сферах деятельности. При этом все более проявляется связь между достижениями страны в области создания и овладения передовыми компьютерными технологиями и ее возможностями решать ключевые проблемы научно-технического прогресса.
Создание и применение современных супер-ЭВМ стало одним из основных направлений в ведущих странах мира. Производительность суперЭВМ на протяжении последних 20-30 лет возрастала ориентировочно на порядок за каждое пятилетие. Это отчетливо прослеживается по всей совокупности выпускаемой продукции, и нет оснований сомневаться в продолжении этой закономерности на следующие годы. К апрелю 2010 г. ИПМ имени М.В. Келдыша и НИИ «Квант» создали принципиально новый тип гибридного вычислительного комплекса и одноименный суперкомпьютер МВС-Экспресс [12], состоящий из 8 узлов, объединенных полносвязным коммутатором PCI-Express. В состав каждого вычислительного узла входят 2 четырехъядерных процессора Opteron с частотой 2.6 ГГц и оперативной памятью 16 ГБ. Производительность вычислительного комплекса составляет 6 терафлоп.
Задача эффективной реализации параллельных алгоритмов всегда считалась непростой хотя бы потому, что программная реализация должна учитывать архитектурные особенности конкретной вычислительной системы. Архитектура вычислительной системы МВС-Экспресс представлена на рис.3.4. суперкомпьютер К-100 состоит из 64 узлов-модулей. В каждом из них по два шестиядерных процессора Хеоп Х5670 фирмы Intel и три графических ускорителя NVIDIA Fermi С2050. В каждом из этих ускорителей 400 арифметических устройств, есть собственная память. Общий объем оперативной памяти составляет 16 терабайт. Структурная схема гибридного вычислительного кластера К-100 представлена нарис. 3.5.
Критически важна для суперкомпьютеров скорость передачи информации между процессорами. В новой машине, помимо общепринятой системы передачи данных Infiniband, используется система собственной разработки МВС-Экспресс, которая обладает преимуществами при передаче относительно коротких фрагментов данных.
Данная ЭВМ потребляет очень мало электроэнергии. Если "Ломоносов" расходует мегаватты, то К-100 потребляет около 80 киловатт. Ее суммарная стоимость вместе со вспомогательным оборудованием составила 65 млн. рублей.
В данной работе рассматриваются задачи моделирования физических-систем методами молекулярной динамики. Существенное развитие в этой области связано с использованием параллельных расчетов: либо одной большой системы (по частям), либо множества маленьких систем (одновременно).
Большая часть времени в МД-моделировании тратится на расчет парных взаимодействий: число пар частиц равно N(N — 1)/2, тогда как число операций для интегрирования уравнений движения пропорционально N. Поэтому для достижения высокой производительности достаточно реализовать на GPU параллельный расчет сил взаимодействий.
В работе разработаны и реализованы алгоритмы параллельных расчетов на гибридных вычислительных комплексах при моделировании установления термодинамического равновесия систем веществ в разных агрегатных состояниях. Параллельная реализация адаптирована для применения технологии CUDA с учетом эффективного распределения данных.
Наиболее вычислительно сложным этапом задачи является подсчет сил взаимодействия и потенциальных энергий частиц системы. В основе параллельного алгоритма лежит модель «ящиков», описанная в пункте 3.1.
Параллелизм на данном этапе состоит в том, что разные «ящики» можно обрабатывать параллельно и независимо друг от друга. Чем больше число обрабатываемых «ящиков», тем полнее получается загрузить GPU работой.
Разработанный параллельный алгоритм формирует сетку блоков для GPU следующим образом. Запускается NUMBOX_X х NUMBOX_Y блоков по NUMBOX_Z нитей в каждом, где NUMBOX_X, NUMBOX_Y, NUMBOX_Z - это количество «ящиков» по осям X, Y и Z соответственно. Таким образом, одна нить обрабатывает один «ящик». Перед работой ядра GPU необходимо скопировать координаты частиц и информацию о распределении частиц по «ящикам» в память GPU. По завершении работы полученные силы взаимодействия и потенциальные энергии частиц копируются из памяти GPU в память CPU. Результаты тестирования параллельного алгоритма приведены в таблице 3.1.
