CC BY
19
УДК 004.942
Об ускорении вычислений при помощи быстрой оценки близости частиц в методе молекулярной динамики
А.А. Харитонов, И.Ф. Ясинский ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина»,
г. Иваново, Российская Федерация E-mail: igor2266@yandex.ru
Авторское резюме
Состояние вопроса: Метод молекулярной динамики активно используется при исследовании процессов в сплошных средах. Его применяют при проектировании ядерных реакторов для изучения процессов воздействия излучения на вещество, в расчетах гидравлического, фильтрационного и температурного режимов для обоснования надежности и безопасности гидротехнических сооружений и при решении ряда других задач энергетики, решаемых с помощью механики сплошных сред. Важным вопросом в методе молекулярной динамики остается выбор способа оценки близости частиц, составляющих моделируемую систему. Критерии выбора такого способа зависят от опыта исследователя, в то время как оценка близости является вычислительно сложной процедурой и значительно сдерживает возможности по изучению больших объектов. Целью исследования является ускорение вычислений в методе молекулярной динамики путем разработки быстрого метода оценки близости частиц. Материалы и методы: Применяются методы молекулярной динамики, способы оценки близости частиц по Эйлеру, Лагранжу. Используются такие методы прикладной математики, как численное интегрирование дифференциальных уравнений методом Рунге-Кутта, оценка погрешности результатов и компьютерное моделирование физических процессов.
Результаты: Предложен комбинированный способ оценки близости частиц, основанный на сортировке частиц по одной из пространственных координат и последовательном отсечении удаленных частиц в соответствии с функцией-метрикой. Выполнено исследование предложенного способа оценки близости частиц на примере физической задачи о динамике газопылевого облака. Численные эксперименты с последовательными и параллельными программными моделями показали увеличение скорости вычислений предложенного подхода по сравнению с традиционными схемами близости по Эйлеру и Лагранжу от 2,2 раз в системе из 103 частиц и до 200 раз в системе из 106 частиц.
Выводы: Предложенный подход при оценке близости частиц в молекулярной динамике позволяет уменьшить временные затраты при моделировании молекулярных процессов и ускорить исследование больших систем, состоящих из более 106 частиц. Способ является альтернативой известным методам Эйлера и Лагранжа, совместим со всеми видами потенциальных взаимодействий и может быть использован при математическом моделировании на основе метода молекулярной динамики. Достоверность результатов подтверждается подобием тестовых и изучаемых в энергетике моделей.
Ключевые слова: методы молекулярной динамики, оценка близости, параллельные вычисления, моделирование молекулярных процессов, алгоритм сортировки.
On accelerating computation through quick evaluation of particle proximity in the method of molecular dynamics
A.A. Kharitonov, I.F. Yasinsky Ivanovo State Power Engineering University, Ivanovo, Russian Federation E-mail: igor2266@yandex.ru
Abstract
Background: The method of molecular dynamics (MMD) is widely applied to studying continuous media processes. It is used in the design of nuclear reactors to analyse the effects of radiation on matter, in the calculation of hydraulic, filtration and temperature conditions to justify the reliability and safety of hydraulic structures, and in a number of other energy problems solved by the mechanics of continuous media. An important MMD problem that has to be addressed is the choice of a method for estimating the proximity of the particles making up the simulated system. The selection criteria for this method depend on the researcher's experience, while proximity estimation is a computationally complex procedure, and it significantly inhibits the ability to study large objects. The aim of the study is to accelerate calculations in the molecular dynamics method by developing an improved technique for estimating the particle proximity. Materials and methods: We used methods of molecular dynamics, Eulerian and Lagrangian particle proximity evalution methods as well as methods of applied mathematics such as numerical integration of differential equations by the Runge-Kutta method, estimation of the error in the results and computer modeling of physical processes. Results: A combined method for estimating the proximity of particles is proposed based on the sorting of particles along one of the spatial coordinates and sequential clipping of the removed particles in accordance with the metric function. The proposed method for estimating the proximity of particles is studied using the example of the physical problem of gas-dust cloud dynamics. Numerical experiments with sequential and parallel software models have shown an speed increase in the computations made in accordance with the proposed approach in comparison with the traditional Eulerian and Lagrangian proximity schemes by 2,2 times in a system of 103 particles and up to 200 times in a system of 106 particles.
Conclusion: A promising approach is proposed for estimating the proximity of particles in molecular dynamics. Reducing the time spent on modeling molecular processes speeds up the study of large systems consisting of more than 106 particles. The method is an alternative to the known Eulerian and Lagrangian methods and is compatible with all kinds of potential interactions. It can be used in many applied fields, in particular in the molecular dynamics mathematical modeling. The reliability of the results is confirmed by the similarity of the test models and the models studied in power engineering.
Key words: molecular dynamics, proximity estimation, parallel computations, molecular process simulation, sorting algorithm.
DOI: 10.17588/2072-2672.2017.4.050-055
Введение. Методы частиц и молекулярной динамики активно используются в данное время в различных областях науки. Нарушение континуальности материалов при сильном деформировании и разрушении создает серьезные сложности в описании подобных процессов в рамках классической механики сплошной среды. Практически любые наноструктуры могут быть смоделированы с чрезвычайно высокой степенью точности на современных многопроцессорных вычислительных системах. Поэтому метод молекулярной динамики (ММД) является важнейшим теоретическим инструментом для разработки нанотехнологий в механике материалов [1].
Для обоснования надежности и безопасности гидротехнических сооружений должны выполняться расчеты гидравлического, фильтрационного и температурного режимов, а также напряженно-деформированного состояния системы «сооружение-основание» на основе применения современных, главным образом, численных методов механики сплошной среды с учетом реальных свойств материалов и пород оснований1.
Описание взаимодействий тяжелых ионов (с энергией до нескольких кэВ) с конденсированными средами основано на ММД. Данный метод позволяет получить значительно больше информации об исследуемой системе (температура, давление, изменение структуры) по сравнению с другими моделями [2]. Изучение воздействия радиации на ГПУ-цирконий выполняется при помощи метода молекулярной динамики [3]. В ряде других задач, актуальных в энергетике, используется ММД [4, 5].
С помощью ММД изучается движение системы частиц, связанных различными видами взаимодействий. В зависимости от решаемых задач, под частицами могут пониматься молекулы или молекулярные фрагменты, отдельные атомы или их группы [6]. Процессы в молекулярных системах описываются при помощи интегрирования дифференциальных уравнений движения классической ньютоновской механики (1):
б2Г "
m,
dt2
(1)
j=1 j
Гидротехнические сооружения ГЭС и ГАЭС. Условия создания. Нормы и требования. Стандарт Организации ОАО РАО «ЕЭС России» СТО17330282.27.140.002-2008, П. 3.1
где /77, - масса /-и частицы; г, - радиус-вектор /-и частицы; Е] - сила, действующая на /'-ю частицу со стороны 7-й.
При использовании методов молекулярной динамики даже в простом случае число решаемых уравнений может составлять около десяти тысяч [7].
При математическом моделировании таких процессов требуется выбрать характер силы взаимодействия частиц и определить оптимальный метод численного решения, соответствующий требуемой точности расчетов. Характерным временем, на котором рассматривается движение частиц, является период от 10 секунд (внутримолекулярные процессы) до двух минут (кон-формационные перестройки молекул в твердых телах) [8].
Размер шага по времени при интегрировании дифференциальных уравнений окажет значительное влияние на качество получаемого решения. Малые значения шага могут привести к такому увеличению времени расчета, что вычислительно сложная задача превратится в вычислительно невозможную, даже при распараллеливании программного кода на многопроцессорных вычислительных системах (МВС). Напротив, большие значения шага будут способствовать как снижению точности решения, так и возникновению неустойчивости.
Целью исследования является ускорение вычислений в методе молекулярной динамики путем разработки усовершенствованного метода оценки близости частиц.
Потенциал Леннард-Джонса позволяет моделировать такие сложные эффекты, как пластичность материала, разрушение, температурное изменение свойств. Рассмотрим систему (рис. 1), в которой взаимодействие частиц определяется потенциалом Леннард-Джонса (2):
m,
dt2 j
j *t
(
B
\n 2
(2)
^ Г (*)"
где А, В - коэффициенты, определяемые свойствами частиц; - расстояние между /'-й и у-й частицами; п1, п2 - параметры модели, влияющие на характер действующих сил. Возможны случаи, когда п1 = 6, п2 = 12 или п1 = 3, п2 = 4 и др.
а)
attraction
б)
Рис. 1. Модель системы в молекулярной динамике: а - схема взаимодействия частиц; б - потенциал Леннард-Джонса
Система уравнений движения для /-й частицы в проекции на ось х будет иметь следующий вид:
г \
B
dt m м
j *i
B
i V1+1 i \ (rj) (r,)
П 2 +1
(x, - X),
(3)
■ = v.
dx¡
где х/, - проекции координаты и скорости частицы на ось х соответственно.
Для проекций на оси у и г уравнения запишутся аналогично. В системе из N частиц для описания движения в трехмерном пространстве будем иметь 6N дифференциальных уравнений.
В молекулярной динамике актуален вопрос производительности вычислительной техники. Для сокращения временных затрат будем отсекать от рассмотрения частицы, находящиеся на значительном расстоянии от данной.
Об оценке близости частиц в методе молекулярной динамики. Известны способы оптимизации вычислительных затрат при расчете близости взаимодействующих частиц (рис. 2) [1].
Близость по Лагранжу. При удалении частиц друг от друга на расстояние Я и более сила их взаимодействия стремится к нулю (рис. 1,б). Параметр Я называется радиусом взаимодействия и используется для определения ближайших частиц, которые будут участвовать при расчете результирующей си-
N
лы ^(рис. 2,а). Эмпирически установлено,
1=1
1
что Я = 10Я0, где Я0 - расстояние, при котором силы отталкивания и притяжения между двумя частицами равны, т.е.
ЧП2+1 '
(4)
, ЧЛ1+1 , ЧЛ
(Г1) (г)
Недостатком метода является необходимость перебора всего количества частиц при расчете расстояния до текущей частицы. Однако значительная экономия временных затрат достигается тем, что для большинства частиц считается только расстояние без воздействующих сил.
Близость по Эйлеру. Модельное пространство делится на ячейки равного размера (рис. 2,б). Каждая ячейка имеет свой номер. Все частицы попадают в свои ячейки исходя из текущих координат. Расчет взаимодействия будет происходить для частиц, находящихся в одной ячейке либо в ближайших соседних. Распределение частиц по ячейкам может происходить согласно следующим выражениям:
(5)
х'-=1 ч £
у'=1+(£
где Х/, У/ - номер ячейки, в которую попала /-я частица; х/, у/ - координаты частицы; Я - линейный размер ячейки.
Данный метод отличается значительными затратами памяти при увеличении точности сетки. Присутствуют дополнительные расходы вычислительных ресурсов на проход пустых ячеек, организацию расположения частиц в ячейках и их корректировку.
Затраты времени при оценке близости согласно Лагранжу (6) и Эйлеру (7) могут составлять соответственно:
^ - NN -1) - №; (6)
Ьв - N, (7)
где N - число частиц в системе.
Схема «Сортировка-покоординатное отсечение». Предлагается схема расчета близости частиц (рис. 3), заключающаяся в следующем алгоритме действий:
1. Производится сортировка частиц по координате одной из осей. Например, при сортировке по х получаем: 6, 3, 5, 7, 1, 2, 4.
2. Рассматриваются частицы, стоящие в сортированном списке от текущей частицы на расстоянии взаимодействия, меньшем R по данной координате. Считаем текущей 5-ю частицу. Тогда под радиус взаимодействия по координате х попадут частицы 6, 3, 7, 1, 2.
3. Оставшиеся частицы также отсекаются по другим осям с учетом расстояния R. Под радиус R по координате у попадают частицы 3, 7, 2.
4. В оставшемся списке отсекаются частицы, расположенные дальше расстояния R, по формуле Евклида с учетом всех осей
Х1 )2 +(У2 - У1)2
Так, из списка ис-
ключается частица 7.
5. Для оставшихся частиц считается взаимодействие. Таким образом, текущая частица 5 будет взаимодействовать с частицами 3, 2.
Рис. 3. Комбинированная схема оценки близости частиц
Численные эксперименты по ускорению оценки близости частиц. Поведение схемы «Сортировка-покоординатное отсечение» было исследовано на задаче моделиро-
вания процессов в газопылевом облаке частиц, которые находятся в силовом поле массивного центрального тела, расположенного в центре системы [9].
Начальными условиями в данной задаче являются следующие характеристики системы: N - число частиц; do - среднее расстояние мет * N
жду частицами; m - масса частицы; -
отношение массы системы к массе массивного Уху
тела; — отношение рандомной составляющей скорости к линейной скорости вращения. Без учета центрального тела на частицу действуют следующие силы: Г = Fgr + Fd/s + ГгеаЫ; (8)
Гдг = -в
т
г^э= -в*
тУУ'. 1Г4 '
ГгеаС = в *
2 11 та
.13
(9) (10) (11)
где Гдг - гравитационная составляющая; Fd/s - диссипативная составляющая; ГгеаС -реактивная составляющая; в - гравитационная постоянная; т - масса частицы; а - равновесное расстояние; г - расстояние между частицами; р - коэффициент диссипации.
Диссипативная составляющая появляется из-за предположения, что нагретая частица интенсивно испускает газообразное вещество. Таким образом, каждая частица окружена газовой оболочкой. Потери энергии при взаимодействии этих оболочек и обусловливают появление диссипативной силы.
При внесении массивного тела массы М возникает центральное поле и на частицы начинает действовать сила:
г=в * тМ, (12)
И2
2
2
где т - масса частицы; Я - расстояние до частицы.
Были проведены серии численных экспериментов с рассмотренной моделью, позволившие оценить эффективность различных схем расчета близости (см. таблицу, рис. 4). Использовался компьютер с двухъядерным процессором, тактовой частотой 2,4 ГГц и технология распараллеливания ОрепМР. В каждом эксперименты выполнено 100 итераций по времени, величина каждой итерации т = 0,1.
Обсуждение. Описанный способ оценки близости частиц позволяет сократить затраты машинного времени при моделировании молекулярных процессов в 2,2 раза в системе из 103 частиц и до 200 раз в системе из 106 частиц и более раз при дальнейшем увеличении числа объектов. Объемы памяти на хранение отсортированного списка невелики. При распараллеливании данного алгоритма вычисли-
тельно сложной процедурой является сортировка. При этом на ее скорость будет значительно влиять число частиц в системе. Наибольшее время на сортировку требуется на первом шаге расчета для инициализации начальных параметров.
По мере работы модели состояние отсортированного массива будет меняться локально и затраты на сортировку окажутся меньше, по сравнению с начальными. Вычислительные затраты данного способа можно оптимизировать за счет выбора быстрого алгоритма сортировки. Считаем, что данный способ более эффективен при работе в интерфейсе МВС с общей памятью, в то время как для ЭВМ-технологии распределенной памяти (MPI) предпочтителен сеточный метод оценки близости.
Результаты численных экспериментов с различными методами оценки близости
Количество Затраты машинного времени, мс
частиц Комбинированная схема по Лагранжу по Эйлеру
Последова- Параллельный Последова- Параллельный Последова- Параллельный
тельный ва- вариант тельный ва- вариант тельный ва- вариант
риант риант риант
1000 12 7 2355 1872 31 16
10000 149 84 237512 160088 858 468
20000 297 250 947347 654924 2976 1435
50000 741 445 ~1,6 часа ~1,4 часа 16395 8159
100000 2229 1386 ~ 6,6 часа ~ 6 часов 64554 30732
200000 5179 4141 == 24,1 часа 248980 122788
500000 12977 7704 ~ 26 минут ~ 12,5 минуты
1000000 26362 18543 ~1,9 часа ~ 0,95 часа
2000000 51945 29133 ~ 8,8 часа ~ 5 часов
10" 10' 10 Рис. 4. Временные затраты при моделировании молекулярной динамики с различными способами оценки близости частиц
Предложенный подход в расчете близости частиц будет интересен всем специалистам, осуществляющим вычислительные эксперименты с моделями на основе молекулярной динамики.
Список литературы
1. Кривцов А.М., Кривцова Н.В. Метод частиц и его использование в механике деформируемого тела // Дальневосточный математический журнал ДВО РАН. -2002. - Т. 3, № 2. - С. 254-276.
2. Непрерывно-атомистическое моделирование процессов взаимодействия тяжелых ионов с металлами на высокопроизводительных вычислительных системах / С.Н. Димова, И.В. Пузынин, Т.П. Пузынина и др. // Труды VII Междунар. конф. «Распределенные вычисления и грид-технологии в науке и образовании» (ГРИД 2016). - Дубна, Россия, 2017. CEUR-WS.org
3. Капустин П.Е. Моделирование ГПУ-циркония методом молекулярной динамики // Известия Самарского научного центра Российской Академии Наук. - 2013. -Т. 15, № 4-5. - С. 1131-1136.
4. Методы молекулярной динамики для моделирования физических и биологических процессов / Х.Т. Хол-муродов, В.М. Алтайский, И.В. Пузынин и др. // Физика элементарных частиц и атомного ядра. - 2003. - Т. 34, вып. 2.
5. Rahman A. Correlations in the motion of atoms in liquid argon // Phys. Rev. - 1964. - 136A. - Р. 405-411.
6. Haile J.M. Molecular dynamics simulation. - Wiley,
1992.
7. Шайтан К. В. Молекулярная динамика пептидов // Динамические модели процессов в клетках и субклеточных наноструктурах / под ред. Г.Ю. Ризниченко и А.Б. Рубина; НИЦ «Регуляция и хаотическая динамика»; Ин-т компьютерных исследований. - М.; Ижевск, 2009. - С. 101.
8. Крокетон К.А. Физика жидкого состояния: пер. с англ. - М., 1978.
9. Динамическая модель образования системы Земля-Луна / Э.М. Галимов, А.М. Кривцов, А.В. Забродин и др. // Геохимия. - 2005. - № 11. - С. 1137-1149.
References
1. Krivtsov, A.M., Krivtsova, N.V. Metod chastits i ego ispol'zovanie v mekhanike deformiruemogo tela [Particle method and its use in deformable body mechanics]. Dal'nevos-tochnyy matematicheskiy zhurnal DVO RAN, 2002, vol. 3, no. 2, pp. 254-276.
2. Dimova, S., Puzynin, I., Puzynina, T., Tukhliev, Z., Hristov, I., Hristova, R., Chernogorova, T., Sharipov, Z. Nepre-ryvno-atomisticheskoe modelirovanie protsessov vzaimodeyst-viya tyazhelykh ionov s metallami na vysokoproizvoditel'nykh vychislitel'nykh sistemakh [Continuum-atomistic modeling of interaction processes of heavy ions with metals by using high performance computer systems]. Trudy VII Mezhdunarodnoy konferentsii «Raspredelennye vychisleniya i grid-tekhnologii v nauke i obrazovanii» (GRID 2016) [Proceedings of the 7th International Conference «Distributed Computing and Grid-technologies in Science and Education» (GRID 2016)]. Dubna, Russia, 2017. CEUR Workshop Proceedings (CEUR-WS.org).
3. Kapustin, P.E. Modelirovanie GPU-tsirkoniya meto-dom molekulyarnoy dinamiki [Modeling of zinconium HCP by the molecular dynamics method]. Izvestiya Samarskogo nauchnogo tsentra Rossiyskoy Akademii Nauk, 2013, vol. 15, no. 4-5, pp. 1131-1136.
4. Kholmurodov, Kh.T., Altayskiy, V.M., Puzynin, I.V., Dardin, T., Filatov, F.P. Metody molekulyarnoy dinamiki dlya modelirovaniya fizicheskikh i biologicheskikh protsessov [Methods of molecular dynamics for modeling physical and biological processes]. Fizika elementarnykh chastits i atomnogo yadra, 2003, vol. 34, issue 2.
5. Rahman, A. Correlations in the motion of atoms in liquid argon. Phys. Rev., 1964, 136A, pp. 405-411.
6. Haile, J.M. Molecular dynamics simulation. Wiley,
1992.
7. Shaytan, K.V. Molekulyarnaya dinamika peptidov [Molecular dynamics of peptides]. Dinamicheskie modeli prot-sessov v kletkakh i subkletochnykh nanostrukturakh [Dynamic models of processes in cells and subcellular nanostructures]. Moscow; Izhevsk, 2009, p. 101.
8. Kroketon, K.A. Fizika zhidkogo sostoyaniya [Physics of liquid state]. Moscow, 1978.
9. Galimov, E.M., Krivtsov, A.M., Zabrodin, A.V., Leg-kostupov, M.S., Eneev, T.M., Sidorov, Y.I. Dinamicheskaya model' obrazovaniya sistemy Zemlya-Luna [A dynamic model of the Earth-Moon system formation]. Geokhimiya, 2005, no. 11, pp. 1137-1149.
Ясинский Игорь Федорович,
ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина»,
кандидат технических наук, доцент кафедры высокопроизводительных вычислительных систем,
e-mail: igor2266@yandex.ru
Yasinsky Igor Fedorovich,
Ivanovo State Power Engineering University,
Candidate of Engineering Sciences (PhD), Associate Professor of the Department of High-Performance Computer Systems, e-mail: igor2266@yandex.ru
Харитонов Александр Андреевич,
ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина»,
магистрант кафедры высокопроизводительных вычислительных систем,
e-mail: b9310508@yandex.ru
Kharitonov Aleksandr Andreyevich,
Ivanovo State Power Engineering University,
Master Course Student of the Department of Department of High-Performance Computer Systems, e-mail: b9310508@yandex.ru
