автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Моделирование теплофизических свойств веществ методами молекулярной динамики с использованием параллельных вычислений

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М.В. Ломоносова

Факультет вычислительной математики и кибернетики

На правах рукописи УДК 519.6

Подрыга Виктория Олеговна

Моделирование теплофизических свойств веществ методами молекулярной динамики с использованием параллельных вычислений

Специальность 05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

1 О НОЯ 2011

Москва-2011

4859201

4859281

Работа выполнена в Институте прикладной математики имени М. В. Келдыша Российской академии наук.

Научный руководитель:

доктор физико-математических наук, профессор,

член-корреспондент РАН Четверушкин Борис Николаевич

Официальные оппоненты:

доктор физико-математических наук, профессор

Попов Александр Михайлович

доктор физико-математических наук Поляков Сергей Владимирович

Ведущая организация:

Московский физико-технический институт (государственный университет).

Защита диссертации состоится » 2011 г. в 15 час. 30 мин. на

заседании диссертационного совета Д 501.001.43 при Московском государственном университете имени М.В.Ломоносова по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет ВМК, аудитория 685.

С диссертацией можно ознакомиться в библиотеке факультета ВМК МГУ имени М. В. Ломоносова. С текстом автореферата можно ознакомиться на официальном сайте ВМК МГУ имени М. В. Ломоносова http://www.cmc.msu.ru в разделе «Наука» - «Работа диссертационных советов» - «Д 501.001.43».

Автореферат разослан « » 2011 г.

Ученый секретарь

диссертационного совета Д 501.001.43 доктор физико-математических наук, профессор

Захаров Е. В.

Общая характеристика работы

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

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

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

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

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

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

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

Цели работы.

1. Моделирование процесса установления термодинамического равновесия систем веществ, имеющих разные структуры, методами молекулярной динамики.

2. Исследование зависимости структуры и характеристик рассматриваемого вещества от геометрии расчетного объема, количества частиц и температуры.

3. Разработка алгоритмов моделирования и анализа результатов молекулярных расчетов для систем атомов в разных агрегатных состояниях.

4. Создание программной реализации с использованием высокопроизводительных технологий.

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

Теоретическая и практическая значимость.

1. Создана оригинальная молекулярно-динамическая программа, позволяющая изучать зависимость структуры и характеристик исследуемого вещества от геометрии расчетного объема, разного количества частиц и от температуры.

2. Получены термодинамические характеристики для систем аргонового газа в объемной геометрии и алюминия в геометрии пластины.

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

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

5. Проведен анализ временного параметра алгоритма теплового контроля (термостата) для эффективного получения желаемой температуры в образце.

6. Разработаны и реализованы алгоритмы параллельных расчетов на гибридных вычислительных комплексах при моделировании установления термодинамического равновесия систем веществ в разных агрегатных состояниях. Параллельная реализация адаптирована для применения технологии США с учетом эффективного распределения данных.

Апробация работы. Основные результаты диссертационной работы докладывались и обсуждались на следующих конференции и семинарах:

1. ХУПГ-й Международной научной конференции студентов, аспирантов и молодых ученых «Ломоносов - 2011» в секции «Вычислительная математика и кибернетика», МГУ имени М. В. Ломоносова (Москва, апрель 2011);

2. научном семинаре «Математическое моделирование» под руководством профессора В. Ф. Тишкина и профессора А. А. Кулешова в ИПМ имени М. В. Келдыша РАН (Москва, май 2011);

3. научном семинаре кафедры вычислительных методов факультета ВМК МГУ имени М. В. Ломоносова под руководством профессора А. В. Гулина (Москва, май 2011);

4. научном семинаре кафедры информатики МФТИ под руководством профессора И. Б. Петрова (Москва, сентябрь 2011).

Публикации. Основные результаты диссертации опубликованы в трех работах, две из которых в изданиях, рекомендованных ВАК [1, 3].

Структура и объем диссертации. Диссертация состоит из оглавления, введения, четырех глав, заключения и списка литературы. Полный объем диссертации составляет 114 страниц, включая 21 иллюстрацию и 2 таблицы, пронумерованные по главам. Список литературы содержит 112 наименований.

Содержание работы

Во введении формулируются цели и задачи диссертации, показывается актуальность, научная новизна, теоретическая и практическая значимость работы. Описывается структура диссертации.

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

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

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

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

В молекулярной динамике частицы двигаются согласно законам классической механики:

т-я=¥ (2.1)

где т — масса частицы, а - ускорение частицы, К - суммарная сила, действующая на данную частицу.

Применив замену вектора ускорения первой производной по времени вектора скорости V, для каждой частицы, значение которого соответствует первой производной по времени от величины радиус-вектора /-ой частицы г() получим систему из 2N обыкновенных дифференциальных уравнений движения для I частиц с массой т, решая которую получим характеристики динамики частиц:

, , / = 1..ЛГ, (2.2)

Л 1

где Р, - равнодействующая сил для /-ой частицы, N - общее число частиц, составляющих систему.

Правая часть уравнения движения (2.1) представляет собой суммарную силу взаимодействия, действующую на частицу с номером /. Силы выражаются через функцию, зависящую от потенциальной энергии.

Потенциальная энергия системы II вычисляется через сумму потенциальной энергии каждой частицы и,:

^ = Ь](и1) (2.3)

1=1

Кинетическая энергия системы КЕ вычисляется как сумма кинетической энергии каждой частицы:

= (2.4)

ш

, «м2

—2~> < = 1"ЛГ.

где ке{ - кинетическая энергия частицы с номером I; |У;| - длина вектора скорости г-ой частицы.

Полная энергия системы Е представляет собой сумму кинетической энергии системы и потенциальной энергии системы:

Е=КЕ+и (2.5)

В первой части главы 2 рассмотрена постановка для системы атомов аргона. Целью задачи является получить равновесное состояние системы

4

частиц с расчетной областью в объемной геометрии, что означает определение

г- 4 г- Го о 5 °о<

< о 4 >°о< о 9. Ь °о о 4 )

< о > °СУ -} ' ^ о ' г. 4 [)

/<1 ^яьс^ / Ч/

7

№4) для £ <х<2£г

(2.6.1)

а б

Рис. 2.1. Графическая реализация периодических граничных условий: а -частицы из рассматриваемой области закрашены в темный цвет [1]; б -распределение атомов алюминия в ГЦК решетке для одной ячейки.

Для моделирования бесконечной области вдоль координатных осей X, У, 2 используются периодические граничные условия с периодами 1Х, Ьу, Ь2 соответственно.

Для частицы, покидающую расчетную область через правую границу х = Ьх:

У = У

х'={х-Ь:с) е

У'=У г' = г

Аналогично для частицы, покидающей расчетную область через левую границу х=0: Гу=У

х'=(х+Ьх) е [0,4) для _4<х<0 (2.6.2)

У'=У

г' = г

Взаимодействие с частицами, находящимися за пределами расчетной области Ьх<х'<(Ьх+гс), моделируется с использованием частиц 0 < х < гс из расчетной области, радиус-векторы которых корректируются следующим образом:

для частиц из расчетной области: г.'=г, ,0<х,<4,

для частиц за пределами расчетной области: ^ = , Ьх<Х1' <{Ьх + гс),

здесь г. - радиус-вектор соответствующей частицы из расчетной области

0 <х{<гс; гс - радиус «обрезания» используемого модельного потенциала.

Аналогичным образом определяются периодические граничные условия вдоль координатных осей У и 2.

Начальные условия включают распределение координат и скоростей частиц.

Рассматривается область кубической формы, заполненная равномерно по осям частицами, количество которых рассчитывается следующим образом:

^Л^.Л'з, (2.7)

где ЛГ], И2, Щ - число частиц по каждой из осей X, У, 2 соответственно.

При условии равномерного заполнения области по осям с одинаковым шагом ДА имеем:

123 Щ ыг Л^з

Тогда общее количество частиц (2.7): Л^л^)3 (2.8)

Координата каждой частицы вычисляется по формуле:

, ] = (2.9)

г°=(0;0;0) , / = 0..(^-1), /Ы1;2;3

Начальные скорости частиц задавались с одинаковыми значениями по абсолютной величине |У0[ и случайно распределенными по направлениям.

Компоненты вектора скорости У(Ух,УуУг) одной частицы получаются следующим образом: Ух =]Уо|-СО80-СО8р,

^=|Уо|.со80-5шр, (2.10)

Гг=|Уо|.5Ш0,

где случайная переменная в равномерно распределена в [-я72; я72), а случайная переменная <р - в интервале [0; 2л).

Выбран парный потенциал взаимодействия, потому межмолекулярный потенциал вычисляется как сумма потенциалов между каждыми двумя частицами с г у =|г; -1\| -расстоянием между /-ой и у'-ой частицами:

и< = £и(гр), У = 1»ЛГ, (2.11)

У

Потенциальная энергия системы С/: ] У

В качестве потенциала взаимодействия выбран потенциал Леннарда-Джонса [2].

V ' \б"

а_

\гч)

- 2

(2.13)

где О - энергия связи молекул (глубина потенциальной ямы в точке минимума потенциальной энергии); а = $2 с - «длина связи» (положение равновесия): ¡и(а)|=£>, к'(а) = 0; а - газокинетический диаметр молекулы. Для аргона эти величины равны: £ = 165*10-23 Дж, а =0,38 нм, ст=0,34 нм.

С целью оптимальности расчетов рассматривается потенциал (2.13) в форме:

«ы=

£

а_

2\ —

г»)

О

"си, . Г..<ГС ^ п

ГЧ>Гс

V2

-2

(2.15)

Здесь и^, способствует сглаживанию функции, то есть помогает избежать разрывность функции при гу=гс; гс = 2,5*ст - радиус «обрезания» потенциала, для аргона гс =0,85 нм.

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

Силы, действующие на частицу с номером /, выражаются через потенциальную энергию этой частицы:

где /(гу) - силы взаимодействия между двумя частицами, которые для потенциала Леннарда-Джонса (2.13) принимают вид:

\13 , ч7"1

121)

(2.17)

Моделирование движения атомов выполнено при помощи разностной схемы Верле [3].

гп,\=гп+у„А(+а"(А1)2 ^ £Л8)

1 _ уп + (а"+1 + а")А? (219)

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

Для решения поставленной задачи используем следующий алгоритм:

1) Расчет по формуле (2.18) новых координат частиц г"+1 для момента времени 1"+1 =г" +Дг, где п - номер шага.

2) Проверка граничных условий (2.6).

3) Вычисление значения парциальной и"+1 (2.11) и суммарной потенциальной энергии системы ип+1 (2.12) частиц и сил (2.16), действующих на частицы.

4) Расчет новых скоростей частиц V"*1 для момента времени /л+1 по выражению (2.19).

5) Вычисление кинетической энергии системы КЕ"'Л по формуле (2.4).

6) Вычисление полной энергии системы Еп+{ (2.5) для момента С*1.

Этапы 1) - 6) повторяются на каждом шаге по времени.

7

С целью слежения за изменением полной энергии системы в работе предлагается алгоритм, в котором вычисляется на каждом шаге величину 8е =\е**1 -ЕпЩЕ'\, характеризующую относительное изменение полной энергии

за время одного шага. Введем ограничение на эту величину сверху и снизу и свяжем ее с изменением временного шага согласно следующему правилу:

1) Если 5£>10"4, то Дг = Дг/2 и пересчитываем все величины с новым шагом.

2) Если 8е <1(Г6, то Дг = ЗД?/2 и пересчитываем все величины с этим новым шагом.

3) Если 10"6 <; 8е < 10"4, то продолжаем счет программы с прежним шагом.

Необходимо отдельно отметить присутствие в знаменателе величины Е".

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

Во второй части главы 2 описана постановка задачи для системы атомов алюминия. Целью задачи является получить равновесное состояние системы при необходимой температуре. Расчетная область имеет геометрию бесконечной пластины, что означает определение по осям X и У периодичных граничных условий с периодами Ьх и Ьу соответственно, по оси 2 условия не накладываются.

Начальное распределение координат определяется в соответствии со структурой кристалла алюминия. Алюминий имеет кубическую гранецентрированную кристаллическую решетку (ГЦК) с параметром (ребром элементарной ячейки) гсг = 0.4 нм.

Рассматриваемая система представляет собой параллелепипед со сторонами Ьх, Ьу, &.

Ь=(1х,Ьу,Ьг)=к-гсг - размеры рассматриваемой области по осям Х,У,2, где к=[кх,ку,к2) - количество элементарных ячеек по направлениям X, 7,1 соответственно.

Таким образом, имеем параллелепипед из кх-ку-кг кристаллов алюминия, решетка кристалла гранецентрированная, что соответствует количеству частиц:

где ^ = кх ку (к2+1) - количество частиц в узлах решетки;

= к* ку (к2+1)+2кх ку кг - количество частиц в центрах граней решетки.

Общее количество частиц: ЛГ=2кх ку(2к2+1) (2.20)

Учитывая граничные условия и используя (2.20), получаем для одного кристалла количество частиц N-6, что можно увидеть на рис. 2.1 б и посчитать по формуле: ¿У=2-М-(2-1+1)=2-3=6.

Координаты задаются следующим образом:

Xi=Xr+l-rcr, ¿ = lJV,/ = 0..(/„-l), л = 1;2;3 (2.22)

где хт1а,ушn,zlmn - точки начала отсчета координат; j\,j2J} - количества точек по направлениям X, Y, Z соответственно.

Значения коэффициентов:

1) для частиц в узлах решетки: = (0, 0, 0), j = (kx,ky,k2 +l), г = 1..Л^, ^,=kxky(kz+l)

2) для частиц в центрах граней, параллельных плоскости (х,у): X^=(6,ö, 0), 8=rcr! 2, j = (kx,ky,kz +1), i = {Nl+\).Jflt N^Nl+kxky{kz +1)

3) для частиц в центрах граней, параллельных плоскости (x.z):

Х** =(0,0,6), Steril, j = (kx,ky,k,), i = (N2+l)..N}, N,=N2 + kxkyk2

4) для частиц в центрах граней, параллельных плоскости(у,г): Х°"=(0,8,5), 5=rcrl2, j = (kx>ky,kz), / = (JV1+I).JV4, N4=N3 +кх ку kz, где N4 =N.

Значения начальных векторов скоростей генерируются из распределения Максвелла, соответствующего значению температуры Г = 10 К. Компоненты вектора скорости V(Vx,VyVz) получены с помощью преобразования Бокса-Мюллера [4]:

Vx = 7-2<t2ln^ cos(27t£),

(2.23)

Fz=A/-2o-2ln^cos(2

где переменные случайные независимые переменные,

равномерно распределенные в (0;1); а1 =кьТ/т - параметр распределения; т -масса частицы, для алюминия т=М-тр =6,8-10"" кг; кь = 1,38-10"23 Дж/К -постоянная Больцмана.

Межмолекулярное взаимодействие описывается с помощью потенциала погруженного атома [5]:

U=£С(лг.) + Уи(/;у), г = 1.JV ,j = 1..jV (2.24)

где и(г^) = н(| гг - Г] |) = и (г) - потенциал парного взаимодействия между i-ой и у-ой частицами; С(щ) - это энергия ¿-го атома за счет воздействия электронов соседних к атомов, и, - суммарная электронная плотность,

создаваемая атомами системы, окружающими ¡'-ую частицу; п(гк) - это

электронная плотность, создаваемая к -ым атомом в месте расположения г-ого атома.

В качестве функции парного взаимодействия и(г^) и функции погружения 1-го атома С(и,) была выбрана форма, предложенная в работе [6]:

<г„) =

1

--Т~а2

iy j

[чц ~Чс)'°\(щ2-afc )6 +a3fli412

где ге-радиус «обрезания» потенциала, для алюминия гс =0.6875 нм.

При большом удалении частиц друг от друга значение энергии взаимодействия таких частиц становится мало. С целью оптимальности расчетов принято рассматривать потенциал (2.24) в заданной форме (2.25) при условиях: и(г>гс) = 0, п(г>гс) = 0.

Константы а, Ь, с для алюминия принимают следующие значения [6]:

a, = 2.9275228176598036; Ь3 = 14.868297626731845;

а2 = 5.1028014804162156; Ь4 = 1.608095393177309;

а3 = 111.37742236893590; с, = 0.58002942432410864;

b, = 8.1106000931637006; с2 = 8.2981185422063639

Ь2 = -334.57493744623503;

Величина силы взаимодействия зависит от потенциальной энергии и имеет следующий вид:

¿(сад-ф^))

--¿г— (2-26>

Для достижения системой необходимой температуры применяется специальный алгоритм - термостат. В работе выбран термостат Берендсена [7]. Изменение кинетической энергии моделируется путем перемасштабирования скоростей атомов молекулярной системы на каждом шаге:

(2.25)

Л= 5__1

(2.27)

где Я - коэффициент пересчета скоростей, Г(/) - текущее значение температуры (в момент времени г), Ть - температура термостата, ть -характерное время взаимодействия с резервуаром.

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

При интегрировании уравнений движения также используется схема Верле: (2.18), (2.19).

Для решения поставленной задачи используем следующий алгоритм:

1) Расчет по формуле (2.18) новых координат частиц г"41 для момента времени - ? + Д/, где п - номер шага.

2) Проверка граничных условий.

3) Вычисление значения потенциальной энергии системы С/л+1 (2.24) частиц и сил Р"+1 (2.26), действующих на частицы.

4) Расчет новых скоростей частиц У"+| для момента времени /и+1 по выражению (2.19).

5) Вычисление коэффициента А"+1 по алгоритму Берендсена (2.27) и масштабирование скоростей.

6) Вычисление кинетической энергии системы КЕ"*'1 по формуле (2.4).

7) Вычисление полной энергии системы Е"+1 (2.5) для момента .

Этапы 1) - 7) повторяются на каждом шаге по времени.

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

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

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

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

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

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

* «V С»". 4 * * © Ф/j, 1м? 5 -J '1» iy-4 - •9 6 f.S ••' > V i 6 !> * iT«-. wr ¿д^.Ф . <5?;

3 с V* t -. i л •••• ...... • •• «'J + ф К Ц-Ч & *•>, «Л е- -' Чф а4 7 „- *

• „- - « t- <- Г " ft 2 «х* ,м „»•; * ¡9 W - v v-tf * ы й o"« t J " ......-■ 9 .....( . Ц'4- л, УЛ 8 tS»

5

ш в

1 в 7

а б

Рис. 3.1. Периодические граничные условия: а - в ячейке 1 моделируемая область, в ячейках 2-9 образы моделируемой области; б - разбиение моделируемой области на кубические «ящики» со стороной, соответствующей радиусу «обрезания» потенциала.

Во второй части главы 3 описывается проблема, заключенная в быстродействии программных комплексов, рассчитанных на решение систем больших размерностей. Решением такой проблемы является использование параллельных вычислений.

В работе разработаны и реализованы алгоритмы параллельных расчетов на гибридных вычислительных комплексах при моделировании установления термодинамического равновесия систем веществ в разных агрегатных состояниях. Параллельная реализация адаптирована для применения технологии CUDA (Compute Unified Device Architecture -унифицированная вычислительная архитектура устройств) с учетом эффективного распределения данных.

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

Наиболее вычислительно сложным этапом задачи является подсчет сил взаимодействия и потенциальных энергий частиц системы. В основе параллельного алгоритма лежит модель «ящиков», описанная в первой части главы 3.

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

Разработанный параллельный алгоритм формирует сетку блоков для GPU следующим образом. Запускается NUMBOX_X х NUMBOX_Y блоков по NUMBOX_Z нитей в каждом, где NUMBOX_X, NUMBOX_Y, NUMBOX Z - это количество «ящиков» по осям X, У, Z соответственно. Таким образом, одна нить обрабатывает один «ящик». Перед работой ядра GPU необходимо скопировать координаты частиц и информацию о распределении частиц по «ящикам» в память GPU. По завершении работы полученные силы взаимодействия и потенциальные энергии частиц копируются из памяти GPU в память CPU. Результаты тестирования параллельного алгоритма приведены в таблице 3.1.

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

Анализ результатов проведения параллельных расчетов на суперкомпьютерах МВС-Экспресс и К-100 подтверждает высокую эффективность вычислительных систем на основе гибридной архитектуры для решения задач молекулярной динамики и перспективность развития

Количество частиц, N Размерность, к:хкухкг Время выполнения последовательного алгоритма, КР,, 10~2 сек. Время выполнения параллельного алгоритма, tm, 10"2 сек. Ускорение,

2 952 6x6x20 12.068 30.223 0.4

4 200 10x10x10 16.169 23.497 0.69

20 200 10x10x50 84.504 40.627 2.08

181 800 30x30x50 795.381 268.71 2.96

901 800 30x30x250 3934.534 988.58 3.98

1 261 800 30x30x350 5497.431 1229.85 4.47

1 441 800 30x30x400 6326.291 1343.16 4.71

3 505 000 50x50x350 15494.465 2604.11 5.95

4 005 000 50x50x400 17681.408 2833.56 6.24

4 020 000 100x100x100 17260.045 2848.19 6.06

8 020 000 100x100x200 34762.332 5172.97 6.72

Таблица 3.1. Результаты тестирования программного комплекса.

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

В первой части главы 4 рассматривается задача для системы частиц аргона. Расчеты заключались в определении времени установления равновесного распределения по скоростям для различного числа частиц и соответствующей температуры.

Расчетная область в форме куба с ребром 1х = Ьу=Ьг- Юнм.

Для N равномерно распределенных в расчетном кубе атомов аргона фиксировались начальные координаты. Скорости у всех частиц задавались одинаковыми по абсолютной величине (|У0| = 5О м/с) и случайно распределенными по направлениям.

Рис. 4.1 Распределение по квадратам скоростей для количества частиц Ы = 103 на момент времени: а - 432 пс; б - 621 пс. Контурная черная линия соответствует линии тренда для данного распределения._

Рис. 4.2. Распределение частиц по квадратам скоростей для Л' = 2 83 на момент времени 432 пс. Контурная черная линия соответствует линии тренда.

Расчеты проводились для количества частиц 103 и 283 со средним значением шага по времени М = 7-1(Г3 пс в течение 60000 и 90000 шагов, общее врем счета составило 432 пс и 621 пс соответственно.

Результаты вычислений для N =103 частиц в случаях 60000 и 90000 шагов представлены на рис. 4.1 в виде распределений по скоростям. Как можно увидеть, система близка к равновесному состоянию уже на момент времени 432 пс (60000 шагов).

Аналогичный результат получается и для другого числа частиц N = 283 в том же объеме, распределение для которого в случае 60000 шагов представлено на рис. 4.2.

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

С \

2-й, 1

(4.1)

у = 1п

И ЛИ2

где пк - количество частиц с квадратами скоростей в интервале

(|НМГ|2+Д|П2);

\У\ - значение скорости по абсолютной величине;

Д|К|2 - интервал квадрата скорости, полученный путем деления максимального значения квадрата скорости на число рассматриваемых в распределении секторов. Максвелловское распределение по квадратам скоростей можно записать в следующем виде

п„(\ V р,| V р +ДIV Р) = С(Г) ехр[-^]^?рД | V р (4.2)

V У

Тогда функция у(\У\2), определенная выше, отличается от показателя экспоненты в этом распределении лишь константой 1п(С(Г)), которая не влияет на наклон отрезка, по которому возможно оценить температуру системы.

Таким образом, при достижении системой равновесного распределения по скоростям строится график зависимости функции у (4.1) от значений квадратов

скоростей \¥\2. Затем с помощью метода наименьших квадратов на участке с наиболее представительной выборкой определяем линейную зависимость вида

у=Р- И2,

в которой коэффициент ¡3 определяет наклон линии. Значение температуры вычисляется по этому коэффициенту:

т т \ .. „.

Т = -йь~Р ^ Расчет температуры системы, определенной таким образом, на момент времени 432 пс (60000-ый шаг) для числа атомов аргона И = 103, дает значение Т = 88 К, а для числа частиц Дг = 283 Т = 63,1 К.

Другой способ вычисления температуры через приходящуюся на одну частицу среднюю кинетическую энергию:

з(М)

дает следующие результаты:

для ЛГ = 103 Г = 87,2К,для N = 2$ Т = 62 К.

Таким образом, температуры, полученные с помощью метода молекулярной динамики путем усреднения параметров движения классической системы частиц (4.4) и с помощью определения коэффициентов наклона по распределению частиц по скоростям (4.3), близки по своим значениям. Откуда следует, что в рассмотренный момент времени большинство частиц термолизовалось, то есть распределение квадратов скоростей частиц близко к максвелловскому. Однако, профиль распределения при количестве частиц N = 283 более похож на профиль максвелловского распределения, что говорит о том, что система при А' = 283 более сильно стремится к термодинамическому равновесию.

Во второй части главы 4 рассматривается задача для системы частиц алюминия.

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

В рассматриваемой области формы параллелепипеда для N частиц алюминия в ГЦК решетке фиксировались начальные координаты. Скорости у всех частиц задавались по распределению Максвелла, соответствующему небольшому значению температуры (т = 10 К).

В ходе моделирования температура системы повышалась за счет использования термостата Берендсена. Целью данной работы являлось получить равновесное состояние системы при температуре т = 300 К.

я С. N N. й 01 С1 " ^ N п й

№8 231? 4375 ЗЭ 121 £13 8123181 Яй

а б

Рис. 4.3. Расчет для количества частиц N = 2952 на момент времени 40 пс при параметре термостата ть = 1 пс: а - распределение частиц по квадратам скоростей (контурная черная линия соответствует линии тренда для данного распределения); б - график изменения температуры в процессе моделирования.

Вычисления проводились для расчетной области 6хбх20 кристаллов алюминия, что соответствует количеству частиц лг = 2952, и 10x10x50 кристаллов, где N = 20200 атомов. Среднее значение шага по времени было Д? = 2-10~3 пс. Изучалось поведение системы при разных временах взаимодействия с резервуаром (разных значениях временного параметра алгоритма Берендсена), рассматривались времена порядка 1 пс и 1 о пс. Расчеты велись в течение 20000 и 40000 шагов, общее время счета составило 40 пс и 80 пс соответственно.

Результаты вычислений для N = 2952 частиц в случае 20000 шагов при параметре термостата ть =1 пс представлены на рис. 4.3 в виде распределения по скоростям и графика изменения температуры в процессе моделирования, при параметре ть -10 пс - на рис. 4.4 соответственно.

а б

Рис. 4.4. Расчет для количества частиц N = 2952 на момент времени 40 пс при параметре термостата т4=10 пс: а - распределение частиц по квадратам скоростей (контурная черная линия соответствует линии тренда для данного распределения); б - график изменения температуры в процессе моделирования.

а б

Рис. 4.5. Расчет для количества частиц N = 2952 на момент времени 80 пс при параметре термостата г4 = 1 пс: а - распределение частиц по квадратам скоростей (контурная черная линия соответствует линии тренда для данного распределения); б - график изменения температуры в процессе моделирования.

Система находится на переходе в равновесное состояние на момент времени 40 пс, но график температуры имеет большую амплитуду относительно желаемого значения. В случае 40000 шагов были сделаны аналогичные расчеты, которые для параметра термостата ть=\ пс представлены на рис. 4.5. На момент времени 80 пс система близка к равновесному состоянию и колебания значений температуры минимальны - желаемый результат получен. Для параметра термостата ть =10 пс на момент времени 80 пс исследуемые графики изменились незначительно.

Вычисления для количества частиц N = 20200 показали, что равновесное состояние системы и достижение небольших амплитуд относительно необходимой температуры достигаются при параметре термостата ть =1 пс уже на момент расчета 40 пс, результаты представлены на рис. 4.6 в виде распределения по скоростям и графика изменения температуры в процессе моделирования.

Рис. 4.6. Расчет для количества частиц N = 20200 на момент времени 40 пс

при параметре термостата ть =1 пс: а - распределение частиц по квадратам скоростей (контурная черная линия соответствует линии тренда для данного распределения); б - график изменения температуры в процессе моделирования.

Аналогично первой части главы 4 с помощью (4.3) вычисляются значения средних температур, основанные на изучении экспоненциального показателя равновесного распределения. Для числа атомов N = 2952 на момент времени 80 пс (40000-ый шаг) с параметром термостата г4 =1 пс получено значение Г = 303,5 К, а для числа частиц N = 20200 с параметром термостата ть =1 пс на момент времени 40 пс (20000-ый шаг) г = 305 К.

Также была высчитана температура через приходящуюся на одну частицу среднюю кинетическую энергию по выражениям (4.4). Для числа атомов N = 2952 получено значение Т = 302 К, для N = 20200: Т = 304 К.

Таким образом, температуры, полученные с помощью метода молекулярной динамики путем усреднения параметров движения классической системы частиц (4.4) и с помощью определения коэффициентов наклона по распределению частиц по скоростям (4.3), близки по своим значениям. Это означает, что в рассмотренный момент времени большинство частиц термолизовалось, то есть распределение квадратов скоростей частиц близко к максвелловскому. Однако,

18

профиль распределения при количестве частиц К = 20200 более похож на профиль максвелловского распределения, что говорит о том, что система при Л^ = 20200 частиц более сильно стремится к термодинамическому равновесию и более быстро, о чем говорит меньшее значение времени счета.

В заключении сформулированы основные результаты диссертационной работы.

Основные результаты работы

1) Исследованы алгоритмы и математические модели для молекулярно-динамических расчетов.

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

3) Проведено исследование зависимости скорости установления равновесного распределения, структуры и характеристик рассматриваемого вещества от геометрии расчетного объема, количества частиц и температуры.

Список цитируемой литературы

1) Allen М.Р. Introduction to Molecular Dynamics Simulation // Computational Soft Matter: From Synthetic Polymers to Proteins. John von Neumann Institute for Computing, Julich, NIC Series. - 2004 - Vol 23.-P. 1-28.

2) Lennard-Jones J. E. Wave functions of many-electron atoms // Proc. Camb. Phil. Soc. - 1931. - Vol. 27. - P. 469.

3) Verlet L. Computer «experiments» on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules // Phys. Rev. -1967.-Vol. 159.-P. 98-103.

4) Box G. E. P., Muller M. E. A Note on the Generation of Random Normal Deviates // The Annals of Mathematical Statistics. - 1958. - Vol. 29, No. 2.-P. 610-611.

5) Daw M. S. Model of metallic cohesion: The embedded-atom method // Phys. Rev. - 1989. - Vol. 39, No. 11. - P. 7441-7452.

6) Zhakhovskii V. V., Inogamov N. A., Petrov Yu. V., Ashitkov S. I., Nishiharaa K. Molecular dynamics simulation of femtosecond ablation and spallation with different interatomic potentials // Applied Surface Science. -2009. - Vol. 1255. - P. 9592-9596.

7) Berendsen H.J.C., Postma J.P.M., W.F. van Gunsteren et al. Molecular dinamics with coupling to an external bath // J. Chem. Phys. - 1984. -Vol. 81.-P. 3684-3690.

Публикации автора по теме диссертации

Подрыта В.О. Моделирование процесса установления термодинамического равновесия методом молекулярной динамики // Математическое моделирование. - 2010. - Т. 22 , No 11 - С 39-48.

Podryga V.O. Molecular Dynamics Method for Simulation of Thermodynamic Equilibrium. // Mathematical Models and Computer Simulations. - 2011. - Vol. 3, No. 3. - P. 381-388. ISSN 2070_0482 Подрыга В.О. Алгоритм моделирования теплофизических свойств металлов методами молекулярной динамики // XVIII Международная научная конференция студентов, аспирантов и молодых ученых «Ломоносов - 2011», секция «Вычислительная математика и кибернетика»: Москва, МГУ имени М. В. Ломоносова, 11-15 апреля 2011 г.: Сб. тезисов. - М.: МАКС Пресс, 2011.-С. 134-136.

Подрыга В.О. Моделирование процесса установления термодинамического равновесия нагретого металла // Математическое моделирование. - 2011. - Т. 23 , No. 9. - С. 3-17.

Подписано в печать: 18.10.2011

Заказ № 5992 Тираж -100 экз. Печать трафаретная. Типография «11-й ФОРМАТ» ИНН 7726330900 115230, Москва, Варшавское ш., 36 (499) 788-78-56 www.autoreferat.ru

Введение.

Глава 1.

Моделирование задач молекулярной динамики.

1.1 Алгоритмы молекулярной динамики.

1.2 Граничные условия.

1.3 Начальные условия.

• Начальные координаты

• Начальные скорости

1.4 Потенциалы взаимодействия.

1.5 Термостаты.

1.6 Численные методы.

Глава 2.

Постановка задачи установления термодинамического равновесия.

2.1 Постановка задачи для аргона.

2.1.1 Граничные условия.

2.1.2 Начальные условия.

2.1.3 Потенциалы.

2.1.4 Численная схема решения.

2.2 Постановка задачи для алюминия.

2.2.1 Граничные условия.

2.2.2 Начальные условия.

2.2.3 Потенциалы.

2.2.4 Алгоритм нагрева (термостат).

2.2.5 Численная схема решения.

Глава 3.

Программная реализация алгоритма и параллельные вычисления.

3.1 Тонкости программной реализации.

3.2 Параллельные вычисления.

Глава 4.

Физическая задача, анализ данных, результаты вычислений.

4.1 Задача для аргона.

4.2 Задача для алюминия.

Современная вычислительная техника дает возможность с помощью компьютерного моделирования решать недоступные для аналитического решения задачи на микроскопическом уровне. Компьютерная молекулярная динамика является одним из наиболее мощных вычислительных методов, эффективно применяемых для моделирования физических, химических и биологических систем [8, 22, 56, 57, 102]. Метод молекулярной динамики (МД) обладает высоким пространственно-временным разрешением и позволяет получить информацию о процессах, происходящих в атомно-молекулярных масштабах и на временах порядка нескольких наносекунд.

В классической молекулярной динамике [30, 89] исследуемая система представляется совокупностью взаимодействующих атомов или молекул, движущихся по законам классической механики.

Развитие вычислительной техники позволяет моделировать динамику молекулярных систем, состоящих из огромного числа частиц, с большим набором параметров и разнообразных условий, имитирующих физический эксперимент. Благодаря развитию и применению компьютерных технологий и графических методов анализа эффективность применяемых методов МД-моделирования неуклонно возрастает. Моделирование реальных физических систем, например, кристаллов или огромных биологических молекул на базе методов МД представляет собой определяющее перспективное направление в ближайшем будущем. В связи с открытием принципиально новых механических и физических свойств у материалов, имеющих структурные элементы ч нанометрового масштаба, чрезвычайно повысился интерес к моделированию материалов на микроскопическом масштабном уровне [25, 32, 33, 94]. Следует также отметить, что в настоящее время все больший интерес у специалистов по математическому моделированию вызывают гетерогенные вычислительные системы, основанные на графических ускорителях; алгоритмы молекулярной динамики являются одними из алгоритмов, эффективно реализующихся на подобных архитектурах.

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

Метод молекулярной динамики получил свое применение в задачах моделирования пластичности и разрушения твердых тел под воздействием высокоскоростного растяжения [23, 24, 26], моделирования изотермического отжига двухкомпонентных кластеров систем металлов [20], в механике деформации твердых тел [25], в моделировании динамики дислокаций [27], в механизмах формирования и разрушения нанорельефа поверхности под действием фемтосекундных лазерных импульсов [1, 14, 59], в динамике биополимеров [52] и многих других актуальных задачах. Необходимым условием для расчета задач моделирования динамических процессов, решаемых с помощью метода молекулярной динамики, является равновесное макроскопическое состояние системы при выбранных тепловых условиях в начальный момент времени. Характеристики состояния термодинамического равновесия системы чаще всего получают экспериментально, что является дорогостоящим, трудоемким и длительным процессом. Программная реализация разработанных и предложенных в данной работе алгоритмов позволяет решать численно задачи релаксации системы, затратив относительно небольшое время по сравнению с реальными экспериментами. В работе рассмотрены и проанализированы модели молекулярных систем в различных агрегатных состояниях (на примере аргона, алюминия), для которых решены задачи релаксации и получены характеристики термодинамического равновесия. Результатами работы алгоритма, наряду с другими характеристиками, являются распределения координат и скоростей атомов вещества в равновесном термодинамическом состоянии. Полученные распределения могут быть использованы в качестве начальных условий для перечисленных выше задач моделирования молекулярно-динамических систем.

Одним из основных этапов решения задач в молекулярной динамике является выбор модели взаимодействия атомов системы, которая позволит получать в поставленных условиях реалистичные данные, наиболее схожие с результатами эксперимента. Взаимодействие частиц описывается с помощью потенциала, который определяется исходя из свойств моделируемого вещества, его агрегатного и термического состояний. Для расчета сил взаимодействия атомов аргона был использован хорошо зарекомендовавший себя потенциал Леннарда-Джонса, позволяющий с минимальными вычислительными затратами получать достаточно точные результаты для систем, характеризующихся парным взаимодействием, в частности, для инертных газов. Для описания свойств металлов парные взаимодействия не обеспечивают качественных результатов [17, 47, 68], так как при расчете сил, действующих в системе, необходимо учитывать вклад коллективизации электронов проводимости. Не смотря на существование большого количества различных видов многочастичных потенциалов, которые используются для решения рассмотренной проблемы, ограничения в подборе параметров и их неспособности воспроизвести полный набор физических характеристик порождают существенные сложности в выборе подходящего потенциала при моделировании процессов в узконаправленных задачах. В диссертационной работе был проведен анализ многочастичных потенциалов, по результатам которого было установлено, что наиболее соответствующим требованиям задачи является потенциал погруженного атома [65-67] в аналитической форме, предложенной в работе [109]. Данная форма оптимизирована для моделирования свойств алюминия как при стандартных условиях, так и в сильно сжатом или наоборот разреженном состояниях.

В литературе [35, 50, 56, 72, 77, 96] предложен ряд методов расчета траекторий движения частиц, анализа их поведения, изучения термодинамических характеристик системы. Однако, как и для каждого быстро развивающегося научного направления имеется потребность в разработке новых алгоритмов, обладающих более совершенными свойствами.

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

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

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

На каждом шаге моделирования решается система обыкновенных дифференциальных уравнений второго порядка, соответствующая второму закону Ньютона и описывающая движение частиц молекулярно-динамической системы. В работе было проведено исследование численных схем, аппроксимирующих уравнения движения. Выбор схемы численного интегрирования обуславливается требованиями устойчивости и минимальности затрат машинного времени при расчете точек фазовой траектории. Для ряда разностных схем таких, как метод Рунге-Кутта, обладающих высоким порядком точности, необходимо проводить неоднократные вычисления правой части уравнений движения на каждом временном шаге, что поглощает основную часть машинного времени. Чем выше порядок точности разностной схемы, тем она сложнее. Это требует, в свою очередь, большой оперативной памяти и одновременно приводит к увеличению числа операций, которые необходимы для перехода к следующему временному шагу, что снижает эффективность. В задачах молекулярной динамики приходится идти на компромисс между точностью и эффективностью, поэтому была выбрана сравнительно простая в реализации и экономичная схема Верле в скоростной форме [106, 107], которая позволяет статистически правильно воспроизводить поведение элементов в течение большого промежутка времени.

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

Одними из основных инструментов моделирования задач МД в России являются, к сожалению, иностранные пакеты программирования, такие как LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [81], HOOMD (Highly Optimized Object Oriented Molecular Dynamics) [58], GROMACS (GROningen MAchine for Chemical Simulations) [76, 105], NAMD (NAnoscale Molecular Dynamics) [64, 80, 91, 93]. По этой причине актуальным является развитие собственных компьютерных разработок - разработок программных комплексов, позволяющих моделировать задачи молекулярной динамики без использования иностранных источников и являющихся более доступными в использовании. Это являлось одной из основных целей данной работы.

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

Диссертация состоит из введения, четырех глав, заключения и списка литературы.

Основные результаты диссертации состоят в следующем.

1) Исследованы алгоритмы и математические модели для молекулярно-динамических расчетов.

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

3) Проведено исследование зависимости скорости установления равновесного распределения, структуры и характеристик рассматриваемого вещества от геометрии расчетного объема, количества частиц и температуры.

Заключение

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

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

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

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

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

В работе разработаны и реализованы алгоритмы параллельных расчетов на гибридных вычислительных комплексах при моделировании установления термодинамического равновесия систем веществ в разных агрегатных состояниях. Параллельная реализация адаптирована для применения технологии С1ЮА с учетом эффективного распределения данных.

Анализ результатов проведения параллельных расчетов на супер-ЭВМ МВС-Экспресс и К-100 подтверждает высокую эффективность вычислительных систем на основе гибридной архитектуры для решения задач молекулярной динамики и перспективность развития данного направления.

1. Ашкрофт Н., Мермин Н. Физика твердого тела. Т. 1. М.: Мир, 1979.

2. Бокий Г.Б. Кристаллохимия. 3-е, переработанное и дополненное издание. -М.: Наука, 1971.

3. Боресков A.B., Харламов A.A. Основы работы с технологией CUDA. -М.: ДМК Пресс, 2010.

4. Борн М. Атомная физика. М.: Мир, 1970 - 492 с.

5. Бурнышев И.Н., Бесогонов В.В. Моделирование поведения водорода в ОЦК-решетке железа при высоких скоростях нагрева // Сборник докладов ФТТ. Институт физики твердого тела и полупроводников. Беларусь, 2005.

6. Бэдсел Ч., Ленгдон А. Физика плазмы и численное моделирование. — М.: Энергоатомиздат, 1989.

7. Валуев A.A., Норман Г.Э., Подлипчук В.Ю. Метод молекулярной динамики: теория и приложения: Математическое моделирование. Физико-химические свойства вещества / Отв. редакторы Самарский A.A., Калиткин H.H. М.: Наука, 1989.

8. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб: БХВ-Петербург, 2004. - 608 с.

9. Воронова Л.И., Тен Э.А. Учет дальнодействия при распределенном молекулярно-динамическом моделировании систем большой размерности // Электронный журнал «ИССЛЕДОВАНО В РОССИИ». 2004. http://zhurnal.ape.relarn.ru/articles/2004/2068.pdf

10. Вычислительные системы МВС-Экспрес, К-100: http://www.kiam.ru/MVS/research/mvsexpress.html http://www.kiam.ru/MVS/resourses/klOO.html

11. Гулд X., Тобочник Я. Компьютерное моделирование в физике. Часть 1. — М.: Мир, 1990.

12. Жаховский В.В., Иногамов H.A., Nishihara К. Новый механизм формирования нанорельефа поверхности, облученной фемтосекундным лазерным импульсом // Письма в ЖЭТФ.- 2008. Т. 87, Вып. 8. - С. 491-496.

13. Жаховский В.В., Нишихара К., Анисимов С.И., Иногамов H.A. Молекулярно-динамическое волн разрежения в средах с фазовыми переходами // Письма в ЖЭВТ. 2000. - Т. 71, Вып. 4. - С. 241-248.

14. Каплан И.Г. Введение в теорию межмолекулярных взаимодействий. — М.: Наука, 1982.

15. Китайгородскй А.И. Молекулярные кристаллы. — М.: Наука, 1971.

16. Киттель Ч. Введение в физику твердого тела. М.: Наука, 1978.

17. Косилов А.Т., Маливанчук A.A., Михайлов Е.А. Молекулярно-динамическое моделирование двухкомпонентных кластеров Cu-Ni, Cu-Pd // Физика твердого тела. 2008. - Т. 50, Вып. 7. - С. 36-40.

18. Кошкин Н., Васильчикова Е. Элементарная физика. Справочник. — М.: Столетие, 1996.

19. Кривцов A.M. Деформирование и разрушение твердых тел с микроструктурой. М.: ФИЗМАТЛИТ, 2007. - 304 с.

20. Кривцов A.M. Описание пластических эффектов при молекулярно-динамическом моделировании откольного разрушения // Физика твердого тела. 2004. - Т. 46, Вып. 6. - С. 1025-1030.

21. Кривцов A.M., Волковец И.Б., Ткачев И.В., Цаплин В.А. Применение метода динамики частиц для описания высокочастотного разрушения твердых тел. Труды всероссийской конференции «Математика, Механика и Информатика 2002 ». М.: ФИЗМАТЛИТ, 2004. - С. 361-377

22. Кривцов А.М., Кривцова H.B. Метод частиц и его использование в механике деформируемого твердого тела // Дальневосточный математический журнал ДВО РАН. 2002. - Т. 3, No. 2. - С. 254-276.

23. Куксин А.Ю., Стегайлов В.В., Янилкин А.В Атомистическое моделирование пластичности и разрушения нанокристаллической меди при высокоскоростном растяжении // Физика твердого тела. — 2008. — Т. 50, Вып. 11.-С. 1984-1990.

24. Лагарьков А.Н., Сергеев В.М. Метод молекулярной динамики в статистической физике // Успехи физических наук. 1978. - Т. 125, Вып. 3. -С. 409-448.

25. Ландау Л. Д., Лифшиц Е.М. Теоретическая физика. М.: Наука, 1976. -Т. 5. Статистическая физика. Часть 1. - 584 с.

26. Метод молекулярной динамики в физической химии / отв. ред. Товбин Ю.К. М.: Наука, 1996.

27. Мирный В., Фрэнер М. Об одной программе моделирования молекулярной динамики газа с элементами распараллеливания алгоритма // Вычислительные технологии. 2001. - Т. 6, No. 3. — С. 32-50.

28. Морозов И.В. Компьютерное моделирование и суперкомпьютерные технологии // Объединенный институт высоких температур РАН. http://www.ihed.ras.ru/norman/studenty

29. Общая физика. Молекулярная физика: План-конспект семинарских занятий / Краснояр. гос. ун-т; сост. О.И.Москвич, О.Ю.Селиверстова. — Красноярск: РИО КрасГУ, 2006.

30. Подмурная O.A. Определение параметров потенциальной функции парного взаимодействия молекул азота и воды // Письма в ЖТФ. 2004. -Т. 30., Вып. 2.-С. 72-75.

31. Поттер Д. Вычислительные методы в физике. М.: Мир, 1975.

32. Розова М.Г., Шпанченко Р.В. Элементы структурной неорганической химии. Учебно-методическое пособие для студентов 1-ого курса. -Московский Государственный Университет им. М.В. Ломоносова, Химический факультет. М.: МГУ, 2001.

33. Самарский A.A., Гулин A.B. Численные методы. М.: Наука. Главная редакция физ.-мат. Литературы, 1989. - 432 с.

34. Сивухин Д.В. Общий курс физики. М.: Наука, 1990. - Т. 2. Термодинамика и молекулярная физика. - 592 с.

35. Солодовников С.Ф. Справочно-методические материалы к курсу «Основы кристаллохимии». — Новосибирский Государственный Университет. Новосибирск, 2006. - 78 с.

36. Суперкомпьютер "Ломоносов": www.t-platforms.ru

37. Турлей Е.В. Молекулярная динамика и диффузия в биомембранах. Диссертация на соискание ученой степени кандидата физико-математических наук. Московский государственный университет им. М.В. Ломоносова. -М.: МГУ, 2006.-119 с.

38. Турлей Е.В., Шайтан К.В., Балабаев Н.К. Динамическая гетерогенность фосфолипидного бислоя и диффузия молекул на границе раздела фаз // М.: Биофизика. 2005. - Т. 50., Вып. 6. - С. 1042-1047.

39. Усольцева Н.В., Ясинский Ф.Н., Соцкий В.В., Костин М.С. Получение стереоизображений в задачах молекулярной динамики // Вестник ИГЭУ. -2009.-Вып. 4.-С. 1-5.

40. Федоров Е.С. Краткое руководство по кристаллографии. Часть 1. -Спб.: Типогр. Ю.Н. Эрлих, 1891. 98 с.

41. Фейнман Р., Лейтон Р., Сэндс М. Фейнмановские лекции по физике. Т. 1. -М.: Мир, 1967.

42. Харламов B.C., Трушин Ю.В., Журкин Е.Е., Лубов М.Н., Пецольдт Й. Исследование методом молекулярной динамики атомов Si и С и кластеров SiC на поверхности кремния // Журнал технической физики. 2008. - Т. 78, Вып. 11.-С. 104-118.

43. Хеерман Д.В. Методы компьютерного эксперимента в компьютерной физике. М.: Наука, Главная редакция физико-математической литературы, 1990.

44. Химическая Энциклопедия. Т.1. ред. Кнунянц И.Л. - М.: Советская энциклопедия. Москва, 1988.

45. Хокни Р., Иствуд Дж. Численное моделирование методом частиц. -М.: Мир, 1987.-640 с.

46. Холмуродов Х.Т., Алтайский М.В., Пузынин И.В., Дардин Т., Филатов Ф.Т. Метод молекулярной динамики для моделировнаия физических и биологических процессов // Физика элементарных частиц и атомного ядра. -2003. Т. 34, Вып. 2. - С. 474-515.

47. Шайтан К.В., Терешкина К.Б. Молекулярная динамика белков и пептидов. Методическое пособие. М.: Юйкос, 2004. - 103 с. http://www.moldyn.ru/library/manual/

48. Шноль Э.Э., Гривцов А.Г. и.др. Метод молекулярной динамики в физической химии. М.: Наука; 1996.

49. Abinit: http://www.abinit.org/

50. Ackland G. J., Thetford R. An improved N-body semi-empirical model for body-centered cubic transition metals // Philosoph. Mag. A. 1987. - Vol. 56. -P.15.

51. Allen M. P., Tildesley A. K. Computer Simulation of Liquids. Oxford: Clarendon Press, 1987. - 244 p.

52. Allen M.P. Introduction to Molecular Dynamics Simulation // Computational Soft Matter: From Synthetic Polymers to Proteins. John von Neumann Institute for Computing, Julich, NIC Series. 2004. - Vol. 23. - P. 1-28.

53. Anderson J.A., Lorenz C.D. Travesset General purpose molecular dynamics simulations fully implemented on graphics processing units // Journal of Computational Physics. 2008. - No. 227. - P. 5342-5359.

54. Anisimov S. I., Zhakhovskii V. V., Inogamov N. A., Nishihara K., Oparin A. M. and Petrov Yu. V. Destruction of a Solid Film under the Action of Ultrashort Laser Pulse//JETP Letters.-2003.-Vol. 77, No. 11.-P. 606-610.

55. Beeman D. Some multistep methods for use in molecular dynamics calculations // Journal of Computational Physics. 1976. - Vol. 20, Issue 2. -P. 130-139.

56. Berendsen H.J.C., Postma J.P.M., W.F. van Gunsteren et al. Molecular dinamics with coupling to an external bath // J. Chem. Phys. 1984. - Vol. 81. -P. 3684-3690.

57. Berendsen H.J.C., van der Spoel D., van Drunen R. GROMACS: a message-passing parallel molecular dynamics implementation // J. Computer Physics Communications. 1995. - Vol. 91. - P. 43-56.

58. Box G. E. P., Muller M. E. A Note on the Generation of Random Normal Deviates // The Annals of Mathematical Statistics. 1958. - Vol. 29, No. 2. -P. 610-611.

59. Brunner R., Kale L., Phillips J. Flexibility and interoperability in a parallel molecular dynamics code // Object Oriented Methods for Interoperable Scientifc and Engineering Computing. SIAM. 1998. - P. 80-89.

60. Daw M. S. Model of metallic cohesion: The embedded-atom method // Phys. Rev. 1989. - Vol. 39, No. 11. - P. 7441-7452.

61. Daw M.S., Baskes M.I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals // Physical Review B. 1984. — Vol. 29, No. 12. - P. 6443-6453.

62. Daw M. S., Foiles S., Baskes M. I. The embedded-atom method: a review of theory and applications // Materials Science Rep. 1993. - Vol. 9. - P. 251-310.

63. Ercolessi F. A molecular dynamics primer. Spring College in Computational Physics, ICTO, Trieste, 1997.

64. Erkoc S. Empirical many-body potential energy functions used in computer simulations of condensed matter properties // Physics Reports. — 1997. Vol. 278, No. 2.-P. 79-105.

65. Finnis M.W., Sinclair J. E. A simple N-body potential for transition metals // Philosoph. Mag.A. 1984. - Vol. 50. - P. 45-55.

66. Foiles S.M., Baskes M. I., Daw M. S., Embedded-atom method functions for the fee metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys // Physical Review . 1986. -Vol. 33.-P. 7983-7991.

67. Frenkel D., Smit B. Understanding molecular simulation: from algorithms to applications. Academic Press, 2002. - 638 p.

68. GROMACS: Fast, Free and Flexible MD. http://www.gromacs.org

69. Guntert P., Mumenthaler C., Wuthrich K. Torsion Angle Dynamics for NMR Structure Calculation with the New Program DYANA // J. Mol. Biol. 1997. -No. 273.-P. 283-298.

70. Hairer E., Lubich C., Wanner G. Geometric numerical integration illustrated by the Stormer/Verlet method // Acta Numerica, Cambridge University Press. — 2003.- Vol. 12.- P. 399-450.

71. Hairer E., Lubich C., Wanner G. Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations. Second Edition. — Springer Series in Computational Mathematics, 2006.

72. LAMMPS: http://lammps.sandia.gov/

73. Landau L.D., Teller E. On the theory of sound dispersion // Physik. Zeits. Sowjetunion. 1936. - Vol. 10. - P. 34-43.

74. Lennard-Jones J. E. Wave functions of many-electron atoms // Proc. Camb. Phil. Soc. 1931. - Vol. 27. - P. 469.

75. Liu X.-Y., Ercolessi F., Adams J.B. Aluminium interatomic potential from density functional theory calculations with improved stacking fault energy // Model. Simul. Mater. Sci. Eng. 2004. - No. 12. - P. 665-670.

76. Maruyama S. Molecular dynamics method for microscale heat transfer// W. J. Minkowycz, E. M. Sparrow (Eds). Advances in Numerical Heat Transfer. 2000. - Vol. 2, Chap. 6. New York: Taylor & Francis. - P. 189—226.

77. Milstein F. Applicability of exponentially attractive and repulsive interatomic potential function in the cubic crystals // J. Appl. Phys. — 1973. Vol. 44, No. 9. -P.3825-3832.

78. Mishin Y. Embedded-atom potential for B2-NiAl / Mishin Y., Mehl M. J., Papaconstantopoulos D: A. // J. Physical Review B. 2002. - Vol. 65, No. 22. -P. 224114-22412.

79. Mishin Y., Farkas D., Mehl M.J., Papaconstantopoulos D.A. Interatomic potentials for monoatomic metals from experimental data and ab initio calculations // Phys. Rev. 1999. - Vol. 59. - P. 3393-3407.

80. Molecular-Dynamics Simulation of Statistical Mechanical Systems. International School of Physics «Enrico Fermi» / Ed. Ciccotti G., Hoover W.G. -Elsevier Science Publishers B.V., 1986.

81. Morse P. M. Diatomic molecules according to the wave mechanics. II. Vibrational levels // Phys. Rev. 1929. - Vol. 34. - P. 57-64.

82. NVIDIA's Compute Unified Device Architecture (CUD A): http://www.nvidia.com/cuda

83. Phillips J. C., Zheng G., Kumar S., Kale L. V. NAMD: Biomolecular Simulation on Thousands of Processors / Proceedings of the IEEE/ACM SC2002 Conference. IEEE Press, 2002. Technical Paper 277.

84. Raczynsky P., Gburski Z. Molecular dynamics and dielectric relaxation of homocysteine layer between graphite walls computer simulation // Rev. Adv. Mater. Sci.-2010.-No. 23.-P. 175-179.

85. Rahman A. Correlations in the motion of atoms in liquid argon // Phys. Rev. -1964.-Vol. 136.-P. 405.

86. Rapaport D.C. The Art of Molecular dynamics Simulation. Second Edition. -Cambridge University Press, 2004.

87. Ruda M., Farkas D., Abriata J. Embedded-atom interatomic potentials for hydrogen in metals and intermetallic alloys // Phys. Rev. 1996. - Vol. 54, No. 14. -P. 9765-9774.

88. Ruhle V. Berendsen and Nose-Hoover thermostats. 2007. http://www.mpip-mainz.mpg.de/~andrienjkyjournalclub/thermostats.pdf

89. Simonelli G., Pasianot R., Savino E. J. // Phys. Rev. 1994. - Vol. 50, No. 2. -P.727-738.

90. Stillinger F.H., Weber T.A. Computer simulation of local order in condensed phases of silicon//Phys. Rev. B. 1985. - Vol. 31. -P. 5262-5271.

91. Stoddard S. D., Ford J. Numerical Experiments on the Stochastic Behavior ofaLennard-Jones Gas System//Phys. Rev. 1973. - Vol. 8. -P. 1504-1512.

92. Sutmann G. Classical molecular dynamics // In: Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms (eds. J. Grotendorst, et al). Julich: NIC. 2002. - Vol. 10. - P. 211-254.

93. Tersoff J. New empirical approach for the structure and energy of covalent system // Phys. Rev. B. 1988. - Vol. 37, No. 12 - P. 6991-6999.

94. Tersoff J. Empirical interatomic potential for silicon with improved elastic properties // Phys. Rev. B. 1988. - Vol. 38, No. 14. - P. 9902-9905.

95. Van Der Spoel D., Lindahl E., Hess B., Groenhof G., Mark A.E., Berendsen HJ. GROMACS: fast, flexible, and free // J. Computer Chemistry. 2005. -Vol. 26, No. 16.-P. 1701-1718.

96. Verlet L. Computer «experiments» on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules // Phys. Rev. 1967. — Vol. 159. - P. 98103.

97. Verlet L. Computer «experiments» on classical fluids. II. Equilibrium correlation functions // Phys. Rev. 1967. - Vol. 165. - P. 201-204.

98. Williams P. L., Mishin Y., Hamilton J. C. An embedded-atom potential for the Cu-Ag system // Modeling in materials science and Engineering. 2006. — No. 14.-P. 817-833.

99. Zhakhovskii V. V., Inogamov N. A., Petrov Yu. V., Ashitkov S. I., Nishiharaa K. Molecular dynamics simulation of femtosecond ablation and spallation with different interatomic potentials // Applied Surface Science. — 2009. -Vol. 1255.-P. 9592-9596.

100. Zope R.R., Mishin Y. Interatomic potentials for atomistic simulations of the Ti-Al system // J. Physical Review B. 2003. - Vol. 68 - P. 024102.

101. Публикации автора по теме диссертации:

102. По дрыгаВ.О. Моделирование процесса установления термодинамического равновесия методом молекулярной динамики // Математическое моделирование. 2010. - Т. 22 , No. 11. - С. 39-48.

103. V. О. Podryga Molecular Dynamics Method for Simulation of Thermodynamic Equilibrium // Mathematical Models and Computer Simulations. 2011. - Vol. 3, No. 3. - P. 381-388. ISSN 2070JM82

104. ПодрыгаВ.О. Моделирование процесса установления термодинамического равновесия нагретого металла // Математическое моделирование. 2011. - Т. 23 , No. 9. - С. 3-17.