CC BY
145
УЧЕНЫЕ ЗАПИСКИ ПЕТРОЗАВОДСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Март, № 2 Физико-математические науки 2015
УДК 004.052.42:519.686
Дмитрий Сергеевич крупянский
аспирант кафедры физики твердого тела физико-технического факультета, Петрозаводский государственный университет (Петрозаводск, Российская Федерация) krupjanski@rambler. ru
Анатолий Дмитриевич фофанов
доктор физико-математических наук, профессор кафедры физики твердого тела физико-технического факультета, Петрозаводский государственный университет (Петрозаводск, Российская Федерация) afofanov@psu.karelia.ru
О влиянии точности арифметических расчетов на результаты молекулярно-динамического эксперимента*
Рассматриваются результаты молекулярно-динамических экспериментов по кристаллизации атомного кластера MgO. В качестве стартовой конфигурации использовалось случайное распределение 500 ионов магния и 500 ионов кислорода по объему куба со стороной 35А. Моделирование проводилось с использованием потенциала межчастичного взаимодействия в форме Борна - Майера -Хиггенса с интегрированием уравнений движения по схеме Верле. При прочих равных условиях расчеты производились с одинарной и двойной точностью. Для обоих случаев представлены результаты исследования динамики изменения структуры моделируемого объекта. На основе полученных данных сделан вывод о критическом влиянии точности арифметических расчетов на течение процесса упорядочения атомной структуры модельного кластера. Дана рекомендация об использовании двойной точности вычислений.
Ключевые слова: точность вычислений, ошибки округления, молекулярная динамика
введение
Метод молекулярной динамики (МД) является одним из наиболее распространенных вычислительных методов, применяемых для моделирования атомных и молекулярных систем. Суть метода МД состоит в численном решении уравнений движения с использованием ЭВМ [8],
[9], [10]:
dr. dp.
m —- = p., —-dt dt
X f (r).
В совокупности с начальными положениями и скоростями частиц эти уравнения представляют собой задачу Коши, для решения которой наиболее часто используется алгоритм Верле [4]:
rn+1 = К + V-At + 0,5 • an-(At )2, vn+1 = vn + 0,5 •( an+1 + a" )-At.
Значение ускорения i-й частицы получается из второго закона Ньютона ai=Fi (t) / m , где m -масса частицы, а F (t) - суммарная сила, действующая на i-ю частицу со стороны остальных. Сила взаимодействия между частицами вычисляется из потенциала межчастичного взаимодействия по формуле fj = -gradU (rtj).
Результатами проведения МД-эксперимента являются координаты и скорости частиц моде-
© Крупянский Д. С., Фофанов А. Д., 2015
лируемого кластера на каждом временном шаге, а также значения потенциальной и кинетической энергии. Особый интерес для исследования структуры материала представляет равновесная конфигурация атомов, получаемая после завершения эксперимента.
При моделировании системы, состоящей из N частиц, для решения каждого уравнения необходимо рассчитать N (N-1) /2 парных взаимодействий, то есть временные затраты на проведение МД-эксперимента растут пропорционально квадрату количества частиц в системе и для достаточно больших систем могут оказаться неприемлемо высокими.
Одним из наиболее доступных способов сокращения времени моделирования является применение графических процессоров (ГП) для параллельного расчета парных взаимодействий [1], [2], [6]. Использование ГП позволяет получать значительный прирост производительности [3] даже на недорогих моделях графических ускорителей. При этом необходимо помнить об архитектурных особенностях этих вычислительных устройств, таких как падение производительности вычислений при использовании двойной точности.
При проведении расчетов на центральном процессоре (ЦП) ощутимой разницы временных затрат при использовании одинарной (single
О влиянии точности арифметических расчетов на результаты молекулярно-динамического эксперимента
73
precision) и двойной (double precision) точности не обнаруживается. Поэтому вопрос об использовании одинарной точности в МД-расчетах оставался неактуальным. Однако расчеты с двойной точностью на ГП проводятся в несколько раз медленнее, чем с одинарной. Авторами [5], [6],
[12], [14] приводятся данные, свидетельствующие о 2-8-кратном падении производительности при выполнении double-операций на современных графических ускорителях.
В настоящей статье приведены результаты исследования о возможности использования арифметики с одинарной точностью для расчетов методом молекулярной динамики.
МЕТОДИКА ИССЛЕДОВАНИЯ
Для оценки степени влияния точности вычислений на результаты молекулярно-динамического расчета был проведен ряд экспериментов с различными условиями и стартовыми конфигурациями. Объектами исследований являлись наноразмерные кластеры, состоящие из: оксида кремния в аморфном и кристаллическом состоянии, оксида магния, оксида кобальта, а также легированные кобальтом ксерогели жидкого стекла. При использовании одинарной точности в некоторых экспериментах наблюдались изменения структуры моделируемых объектов, отсутствующие в аналогичных экспериментах, проведенных с двойной точностью. Наилучшим образом такие изменения видны в экспериментах, моделирующих фазовый переход со сменой полиморфной модификации или с изменением степени упорядоченности (кристаллизация). При моделировании неупорядоченных структур влияние точности вычислений на поведение структуры атомных кластеров не выявлено из-за сложностей численного описания таких структур. Наблюдаемое в экспериментах появление структурных различий связано с ошибками округления [13], [15], являющимися общей проблемой расчетов с плавающей запятой, не связанной с конкретным типом вычислительных устройств [11].
Влияние точности расчетов на результаты моделирования продемонстрировано на примере наиболее показательного эксперимента, в ходе которого происходит кристаллизация модельного кластера. Результаты экспериментов, полученные с разной точностью, сравнивались с помощью методики, описанной в статье [7] (http://mmp.vestnik. susu.ac.ru/pdf/v7n2st4.pdf). Также в настоящей статье приводятся соответствующие кривые потенциальной энергии, отражающие изменения в структуре моделируемого объекта.
ОПИСАНИЕ ЭКСПЕРИМЕНТА
Стартовая конфигурация представляла собой кластер, имеющий форму куба со стороной 35А и состоящий из 1000 случайно распределенных по объему атомов кислорода и магния в равном стехиометрическом соотношении. Единственное
ограничение, накладываемое на взаимное расположение частиц, заключается в том, что расстояние между двумя атомами не должно быть меньше 1,5А. В расчетах использовался потенциал в форме Борна - Майера - Хиггинза:
J 0 +A, exp
Ра
2
U
iJ
iJ
Исследуемый кластер помещался в центр модельного объема в виде куба со стороной 150А, что было в несколько раз больше размеров исходного атомного кластера. Таким образом, поверхность кластера могла принимать энергетически наиболее выгодную форму. Во всех модельных экспериментах радиус обрезания потенциалов был равен 70А, что позволило учесть взаимодействие каждого атома со всеми остальными атомами моделируемой системы. Один временной шаг соответствовал 10-15с модельного времени, всего было рассчитано 100 000 шагов. В ходе эксперимента происходило зарождение центров кристаллизации с дальнейшим ростом новой фазы по всему объему кластера.
Моделирование проводилось с использованием параллельной реализации метода молекулярной динамики для графического процессора NVIDIA [6] (http:Mibrary.ru/item.asp?id=20318865) с одинарной и двойной точностью на следующей конфигурации:
Intel Core i7-3770L 3500 МГц 8096 Гб DDR3-2133 2ГБ GeForce GTX 660, com-
3.0
g++ 4.6.3, nvcc 5.0 cuda runtime, libquadmath, igraph c library 0.7.1
Полученные результаты были проверены с помощью реализации МД для процессора с архитектурой x86 с четверной точностью (библиотека gcc libquadmath) при прочих равных условиях.
Процессор: ОЗУ:
видеокарта: pute capability Компиляторы: Библиотеки:
результаты
В ходе МД-эксперимента, проведенного с двойной точностью, образуется упорядоченный кластер, имеющий правильную решетку типа NaCl (две ГЦК-подрешетки из катионов и анионов, сдвинутых друг относительно друга на половину телесной диагонали элементарной ячейки) (рис. 1а). При проведении расчетов с одинарной точностью полученный в результате эксперимента кластер представляет собой агрегат из мелких кристаллитов, сросшихся друг с другом под разными углами (рис. 1б). В данном случае более низкая точность расчета влечет за собой потерю вклада удаленных частиц при расчете сил взаимодействия, что негативно сказывается на общей упорядоченности структуры модельного кластера.
74
Д. С. Крупянский, А. Д. Фофанов
Рис. 1. Схематичное изображение структуры модельного кластера, полученного при расчетах с двойной (а) и одинарной (б) точностью
На рис. 2 представлен график изменения потенциальной энергии модельного кластера в ходе экспериментов с одинарной и двойной точностью. Первые 3000 шагов модельного времени кривые полностью совпадают. В момент времени t = 3000 на графике заметно незначительное расхождение кривых, после которого они снова практически совпадают.
Рис. 2. Графики изменения потенциальной энергии при расчетах с разной точностью
Значительное расхождение кривых начинается после 20 000 шагов МД-эксперимента. Величина потенциальной энергии при моделировании с одинарной точностью продолжает равномерно убывать вплоть до окончания эксперимента, а более резкое снижение значений этого параметра при расчете с двойной точностью свидетельс-
твует о существенных изменениях в структуре моделируемого объекта - кристаллизации.
К моменту времени t = 20000 структуры кластеров, полученных с различной точностью, практически неразличимы - оба кластера содержат некоторое количество упорядоченных областей (зародыши кристаллизации). Далее, при использовании двойной точности вычислений один из зародышей становится центром кристаллизации всего кластера. Атомы, окружающие эту область, начинают упорядочиваться, происходит фазовый переход. Так как в расчете сил учитываются взаимодействия каждой частицы со всеми остальными, суммарный вклад удаленных частиц является значимой величиной. Поэтому зародыши, не ставшие центром кристаллизации, под действием силы со стороны уже упорядоченных атомов естественным образом встраиваются в растущую структуру. При использовании одинарной точности вклад одной частицы, расположенной на достаточном расстоянии, может оказаться настолько малым, что будет потерян при округлении. Таким образом, частицы, расположенные на противоположных концах кластера, практически не взаимодействуют друг с другом. В этом случае с момента времени t = 20000 каждый зародыш начинает формировать собственный центр кристаллизации. Близкие центры, взаимодействующие друг с другом, срастаются в один кристаллит, а более дальние срастаются «неправильно», то есть разориентированы друг относительно друга.
Графики на рис. 2 отражают качественные различия в течение процесса моделирования, обусловленные исключительно различной точностью выполнения расчетов. Для более детального анализа этих различий было проведено исследование структурных изменений моделируемого кластера в соответствии с методикой, изложенной в [7]. На каждом шаге МД-эксперимента был построен граф g, описывающий структуру модельного кластера. Вершинам графа g соответствуют найденные в кластере магний-кисло-родные октаэдры, ребро соединяет две вершины, если соответствующие им многогранники имеют
Рис. 3. Графики изменения порядка графа g (слева) и его гигантской компоненты (справа)
О влиянии точности арифметических расчетов на результаты молекулярно-динамического эксперимента
75
общее ребро. Затем с помощью библиотеки igraph был рассчитан порядок гигантской компоненты связности этого графа - инвариант, чувствительный к процессам изменения структуры кластера, происходящим при кристаллизации.
Результаты расчетов инвариантов графа g для одинарной и двойной точности приведены на рис. 3. При расчете с двойной точностью кривые порядка графа g и порядка его гигантской компоненты практически совпадают. Это говорит о том, что в большинстве случаев каждая новая вершина графа g присоединяется к гигантской компоненте. На кривых хорошо виден резкий подъем, соответствующий процессу фазового перехода. Рост кристалла в этом случае происходит из одного центра кристаллизации. При расчете с одинарной точностью порядок графа g растет без значительных скачков, при этом порядок гигантской компоненты этого графа осциллирует около некоторого постоянного значения вплоть до момента времени t = 60000. Таким образом, гигантская компонента в этом случае не формируется. За 80 000 шагов модельного времени образуется около десятка связных компонент, имеющих порядок от 1 до 10 вершин, которые за последние 20 000 шагов эксперимента срастаются в две. Рост кристалла происходит из нескольких центров кристаллизации. Образовавшиеся кристаллиты срастаются под различными углами.
По графикам потенциальной энергии, представленным на рис. 2, можно было бы сделать вывод об адекватности данных, полученных с одинарной точностью за первые 15 000-20 000 шагов модельного времени. Однако графики на рис. 3 свидетельствуют о накоплении заметной ошибки уже за первые 10 000-15 000 шагов.
ВЫВОДЫ
Проведено множество молекулярно-динамических экспериментов с одинарной и двойной точностью арифметических расчетов для наноразмерных кластеров на основе оксида кремния, оксида магния и оксида кобальта. Результаты большинства проведенных вычислительных экспериментов качественно соответствуют результатам эксперимента, описанного в настоящей работе. Собранные данные свидетельствуют о критическом влиянии точности выполнения расчетов на конечную структуру моделируемых атомных кластеров.
На примере кристаллизации кластера на основе оксида магния показано, что использование одинарной точности в молекулярно-динамических расчетах приводит к потере вклада удаленных частиц при расчете сил межатомного взаимодействия, что выражается в снижении степени упорядоченности структуры моделируемого объекта. При моделировании неупорядоченных кластеров выявить различия в ходе экспериментов, проводимых с различной точностью, не удалось в силу сложности численного сравнения структуры этих кластеров.
Результаты всех вычислительных экспериментов, проведенных с двойной точностью, соответствуют результатам аналогичных экспериментов, проведенных с четверной точностью с использованием библиотеки gcc libquadmath. Следовательно, двойную точность для данной задачи и схемы интегрирования можно считать достаточной. Использование одинарной точности расчета может качественно изменить ход вычислительного эксперимента и не может быть оправдано приростом производительности, получаемой на таких вычислительных устройствах, как графические процессоры.
* Работа выполнена при поддержке Программы стратегического развития ПетрГУ в рамках реализации комплекса мероприятий по развитию научно-исследовательской деятельности на 2012-2016 гг.
СПИСОК ЛИТЕРАТУРЫ
1. Боярченков А. С., Поташников С. И. Использование графических процессоров и технологии CUDA для задач молекулярной динамики // Вычислительные методы и программирование. 2009. Т. 10. С. 9-23.
2. Боярченков А. С., Поташников С. И. Параллельная молекулярная динамика с суммированием Эвальда и интегрированием на графических процессорах // Вычислительные методы и программирование. 2009. Т. 10. С. 158-175.
3. Галимов М. Р., Биряльцев Е. В. Некоторые технологические аспекты применения высокопроизводительных вычислений на графических процессорах в прикладных программных системах // Вычислительные методы и программирование. 2010. Т. 11. С. 77-93.
4. Гельчинский Б. Р., Мирзоев А. А., Воронцов А. Г. Вычислительные методы микроскопической теории металлических расплавов и нанокластеров. М.: ФИЗМАТЛИТ, 2011. 200 с.
5. Кривов М. А., Казеннов А. М. Сравнение вычислительных возможностей графических ускорителей NVidia при решении различных классов задач // Труды Всероссийской научно-практической конференции «Применение гибридных высокопроизводительных вычислительных систем для решения научных и инженерных задач». Н. Новгород, 2011. С. 18-24;
6. К р у п я н с к и й Д. С., Л о б о в Д. В., О с а у л е н к о Р. Н. Реализация метода молекулярной динамики посредством технологии Nvidia CUDA // Ученые записки Петрозаводского государственного университета. Сер. «Естественные и технические науки». 2013. № 2 (131). С. 84-86.
7. Крупянский Д. С., Фофанов А. Д. Алгоритм поиска точечных подмножеств и его применение для анализа атомной структуры модельных кластеров // Вестник Южно-Уральского государственного университета. Сер. «Математическое моделирование и программирование». Челябинск, 2014. Т. 7. № 2. С. 46-54.
8. Рит М. Наноконструирование в науке и технике: введение в мир нанорасчета / Пер. с англ. Э. М. Эпштейна. М.; Ижевск: НИЦ «Регулярная и хаотическая динамика», 2005. 160 с.
9. X е е р м а н Д. В. Методы компьютерного эксперимента в теоретической физике / Под ред. С. А. Ахманова. М.: Наука, 1990. 176 с.
76
Д. С. Крупянский, А. Д. Фофанов
10. Холмуродов X. Т., Алтайский М. В., П у зынин И. В., Д а р д и н Т., Филатов Ф. П. Методы молекулярной динамики для моделирования физических и биологических процессов // Физика элементарных частиц и атомного ядра. 2003. Т. 34, вып. 2. С. 472-515.
11. Goldberg D. What every computer scientist should know about floating-point arithmetic // ACM Computing Surveys. 1991. Vol. 23. P. 5-48.
12. Itu L. M., Moldoveanu F., Suciu C., Postelnicu A. Comparison of single and double floating point precision performance for Tesla architecture GPUs // Bulletin of the Transilvania University of Brasov. 2011. Vol. 4 (53). № 2. P. 131-138.
13. NVIDIA Corporation. NVIDIA CUDA C Best Practices Guide ver. 6.0. Santa Clara, Ca, 2014.
14. Wezowicz M., Saunder D., Taufer M. Dealing with performance/portability and performance/accuracy trade-offs in heterogeneous computing systems: A case study with matrix multiplication modulo primes // In Proc. SPIE 8403, Modeling and Simulation for Defense Systems and Applications VII. 2012. P. 8-18.
15. Whitehead N., Fit-Florea A. Precision and Performance: Floating Point and IEEE 754 Compliance for NVIDIA GPUs // Technical white paper by NVIDA. 2011.
Krupyanskiy D. S., Petrozavodsk State University (Petrozavodsk, Russian Federation) Fofanov A. D., Petrozavodsk State University (Petrozavodsk, Russian Federation)
IMPACT OF ARITHMETIC PRECISION ON RESEARCH RESULTS
of molecular-dynamic simulation
Research results of molecular-dynamic simulations, in which crystallization of atomic cluster MgO occurs, are considered in the article. 500 ions of magnesium and 500 ions of oxygen were randomly placed into a cube with a side of 35A as an initial configuration. In the study we used Born - Mayer - Huggins interatomic potential and Verlet scheme for integration of equations of motion. The Simulations were conducted with the single and double precision ceteris paribus. The changes in structural dynamics of the modeled objects are presented for both cases. Based on the obtained results we came to a conclusion that the calculation accuracy significantly influences the progress of ordering the atomic structure of the modeled object. We have shown that molecular-dynamic simulations require strict accuracy provided by double precision.
Key words: calculations accuracy, rounding errors, molecular dynamics
REFERENCES
1. Boyarchenkov A. S., Potashnikov S. I. The Use of GPUs and CUDA Technology for Molecular Dynamics Applications [Ispol’zovanie graficheskikh protsessorov i tekhnologiy CUDA dlya zadach molekulyarnoy dinamiki]. Vychislitel’nye metody iprogrammirovanie. 2009.Vol. 10. P. 9-23.
2. Boyarchenkov A. S., Potashnikov S. I. Parallel Molecular Dinamics with Evald’s Summation and Integration on GPUs [Parallel’naya molekulyarnaya dinamika s summirovaniem Eval’da i integrirovaniem na graficheskikh protsessorakh]. Vychislitel’nye metody iprogrammirovanie. 2009. Vol. 10. P. 158-175.
3. Galimov M. R., Biryal’tsev E. V. Some Technological Aspects of Application of High-performance Calculations on Graphic Processors in Applied Program Systems [Nekotorye tekhnologicheskie aspekty primeneniya vysokoproizvoditel’nykh vychisleniy na graficheskikh protsessorakh v prikladnykh programmnykh sistemakh]. Vychislitel ’nye metody i programmirovanie. 2010. Vol. 11. P. 77-93.
4. Gel’chinskiy B. R., Mirzoev A. A., Vorontsov A. G. Vychislitel ’nye metody mikroskopicheskoy teorii metal-licheskikh rasplavov i nanoklasterov [Computational Methods in Microscopic Theory of Metal Melts and Nanoclasters].
Moscow, 2011. 200 p.
5. Krivov M. A., Kazennov A. M. Comparison of Computing Capabilities of NVidia Graphic Accelerators in the Solution of Various Classes of Tasks [Sravnenie vychisntel’nykh vozmozhnostey graficheskikh uskoriteley NVidia pri reshenii razlichnykh klassov zadach]. Trudy Vserossiyskoy nauchno-prakticheskoy konferentsii "Primenenie gibridnykh vysokoproizvoditel ’nykh vychislitel’nykh sistem dlya resheniya nauchnykh i inzhenernykh zadach” [Proc. All-Russian scientific and practical conference “Application of hybrid high performance computing systems for solving scientific and engineering problems”]. Nizhniy Novgorod, 2011. P. 18-24.
6. Krupyanskiy D. S., Lobov D. V., Osaulenko R. N. Implementation of Molecular Dynamics Method Using Nvidia CUDA Technology [Realizatsiya metoda molekulyarnoy dinamiki posredstvom tekhnologii Nvidia CUDA]. Uchenye zapiski Petrozavodskogo gosudarstvennogo universiteta. Ser. "Estestvennye i tekhnicheskie nauki” [Proceedings of Petrozavodsk State University. Natural and Engineering Sciences]. 2013. № 2 (131). P. 84-86.
7. Krupyanskiy D. S., Fofanov A. D. An Algorithm Searching for Point Subsets with Applications to the Analysis of the Atomic Structure of Modelled Clusters [Algoritm poiska tochechnykh podmnozhestv i ego primenenie dlya analiza atomnoy struktury model’nykh klasterov]. Vestnik Yuzhno-Ural ’skogo gosudarstvennogo universiteta. Ser. "Matematicheskoe modeli-rovanie iprogrammirovanie”. 2014. Vol. 7, № 2. P. 46-54.
8. R i t M. Nanokonstruirovanie v nauke i tekhnike: vvedenie v mir nanorascheta [Nano-Engineering in Science and Technology: An Introduction to the World of Nano-Design]. Moscow; Izhevsk, 2005. 160 p.
9. H e e r m a n D. V. Metody komp’yuternogo eksperimenta v teoreticheskoyfizike [Computer Simulations Method in Theoretical Physics]. Moscow, 1990. 176 p.
10. Kholmurodov Kh. T., Altayskiy M. V., Puzynin I. V., Dardin T., Filatov F. P. Molecular dynamics method for simulations of physical and biological processes [Metody molekulyarnoy dinamiki dlya modelirovaniya fizicheskikh i biologicheskikh protsessov]. Fizika elementarnykh chastits i atomnogoyadra. 2003. Vol. 34, № 2. P. 472-515.
11. Goldberg D. What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys. 1991. Vol. 23. P. 5-48.
12. Itu L. M Moldoveanu F., Suciu C., Postelnicu A. Comparison of single and double floating point precision performance for Tesla architecture GPUs. Bulletin of the Transilvania University of Brasov. 2011. Vol. 4 (53). № 2.P. 131-138.
13. NVIDIA Corp. NVIDIA CUDA C Best Practices Guide ver. 6.0. Santa Clara, CA, 2014.
14. Wezowicz M., Saunder D., Taufer M. Dealing with performance/portability and performance/accuracy trade-offs in heterogeneous computing systems: A case study with matrix multiplication modulo primes. In Proc. SPIE 8403, Modeling and Simulation for Defense Systems and Applications VII. 2012. P. 8-18.
15. Whitehead N., Fit-F lo rea A. Precision and Performance: Floating Point and IEEE 754 Compliance for NVIDIA GPUs. Technical white paper by NVIDA. 2011.
Поступила в редакцию 17.10.2014
