Зачем и какие нужны суперкомпьютеры эксафлопсного класса? Предсказательное моделирование свойств и многомасштабных процессов в материаловедении Текст научной статьи по специальности «Математика»

УДК 519.6

А. Ю. Куксин, А. В. Ланкин, И. В. Морозов, Г. Э. Норман, Н. Д. Орехов, В. В. Писарев, Г. С. Смирнов, С. В. Стариков, В. В. Стегайлов, А. В. Тимофеев

ЗАЧЕМ и КАКИЕ нужны суперкомпьютеры

эксафлопсного класса? Предсказательное моделирование свойств и многомасштабных процессов в материаловедении

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

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

© А. Ю. Куксин, А. В. Ланкин, И. В. Морозов, Г. Э. Норман, Н. Д. Орехов, 2014

© В. В. Писарев, Г. С. Смирнов, С. В. Стариков, В. В. Стегайлов, А. В. Тимофеев, 2014

© Объединенный институт высоких температур РАН, 2014

© Программные системы: теория и приложения, 2014

Введение

В первой части работы рассматривается подход, позволяющий выявить, для каких задач нужны суперкомпьютеры эксафлопсного класса. Имеются в виду задачи, для решения которых требуется загрузка всего или значительной части суперкомпьютера, т.е. сотен тысяч или даже миллионов вычислительных ядер единовременно под одну задачу. Анализ ведётся на примерах из области многомасштабного атомистического моделирования — прорывного направления современной науки. Чтобы стала ясна роль этого направления, укажем, что около трети машинного времени и задач, решаемых на лучших суперкомпьютерах США, принадлежит к нему [1]. При этом основные траты машинного времени приходятся на атомистическое моделирование и квантовые задачи, т.е. на два самых глубоких масштаба. Неудивительно, что, например, при разработке архитектуры суперкомпьютеров IBM BlueGene/L задачи молекулярного моделирования рассматривались в качестве главного приоритета [2].

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

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

(1) Модификация поверхности при облучении металла субпикосе-кундными лазерными импульсами. Используется и развивается двухтемпературная атомистическая модель вещества. Проводятся расчёты, способные напрямую воспроизвести экспериментальные данные по деформации поверхности вследствие процессов плавления и абляции, вызванных лазерным облучением. Развивается полномасштабный трёхмерный подход.

(2) Радиационно-индуцированные структурные изменения в облучённом топливе ядерных реакторов нового поколения на быстрых нейтронах. Развиваются механистические модели: уточняются базовые механизмы диффузии точечных дефектов и их кластеров, диффузии продуктов деления (в первую очередь ксенона), образования кластеров точечных дефектов и пузырей ксенона, рекомбинации точечных дефектов. Определяются микроскопические параметры для указанных моделей.

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

(4) Свойства и структура супрамолекулярных соединений включения— газовых гидратов. Исследуются новые фазы гидратов водорода в области высоких давлений.

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

(6) Теплофизические свойства (теплоёмкость, вязкость и др.) систем пылевых частиц в плазме, влияние анизотропии плазменно-пы-левой системы. Исследуется возможность термодинамического подхода для описания систем малого числа пылевых частиц.

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

(8) Двойной электрический слой на границе углеродного материала и электролита, влияние электронно-дырочной структуры материала электрода на ёмкостные свойства.

(9) Распараллеливание, адаптация программ многомасштабного моделирования для выполнения на гибридных вычислительных

системах, включающих графические ускорители.

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

В заключении даются ответы на вопросы, зачем и какие нужны суперкомпьютеры эксафлопсного класса в научных исследованиях. Их развитие анализируется в [3], где также затронуты требования к архитектуре суперкомпьютеров с точки зрения масштабируемости задач молекулярного моделирования.

1. Зачем нужны суперкомпьютеры эксафлопсного класса?

1.1. Параллелизм классических МД-расчётов

1.1.1. Локальность обменов данными

МД-расчёты проводятся для миллионов и миллиардов частиц [4]. Для распараллеливания расчётный объем разделяется на подобласти (domain decomposition), каждая из которых «поручается» одному ядру, пример одномерного разделения приведён на рис. 1. Ограничение масштабируемости алгоритмов молекулярной динамики определяется межпроцессорным обменом данными, что схематично проиллюстрировано для декомпозиции по пространству в одномерном случае на рис.1 и двумерном — на рис.2. В зависимости от геометрии системы задача разделяется на отдельные пространственные области, объёмы которых считаются на разных процессорах. На границе областей процессоры обмениваются данными для расчёта сил. Пространственная локальность взаимодействия обеспечивает высокую параллельную эффективность данного подхода.

Рис. 1. Схема, поясняющая параллелизацию на основе декомпозиции по пространству в методе МД (короткодействующие потенциалы)

Рис. 2. Схематичное представление локальности обменов данными для двумерной системы при атомистическом моделировании.

Рис. 3. Зависимость ускорения от числа ядер (слева) и оптимальное соотношение числа атомов и ядер(справа)

1.1.2. Масштабируемость

Слева на рис. 3 показана зависимость ускорения от числа ядер для различных размеров систем с потенциалом Леннард-Джонса (указано число частиц в системе). При фиксированном числе частиц увеличение числа подобластей, т.е. увеличение числа используемых вычислительных ядер и уменьшения числа частиц на каждое ядро, вначале даёт линейное ускорение расчётов, а потом проходит через максимум, когда межъядерный обмен становится затратным по времени. Таким образом, для каждого числа вычислительных ядер существует число частиц, оптимальное для расчётов. Справа на рис.3 — оптимальное соотношение числа атомов и ядер. Точками показаны расчёты с эффективностью 80%. Сплошная кривая—ограничение по оперативной памяти. Для сравнения приведён расчёт на суперкомпьютере IBM BlueGene/L [5]. При увеличении числа ядер можно проводить расчёты для большего число частиц и, соответственно, расширяется круг явлений и процессов, доступных для исследования. Такой подход предложен в [5], где он проиллюстрирован на примере пластической деформации и разрушения при высокоскоростном деформировании. Набор других примеров приведён в данной работе. Отметим, что для систем с

механика сплошной среды

квантовые "ab initio" подходы

Рис.4. Ступени многомасштабного подхода и их взаимосвязи

короткодействующими потенциалами типа ЕАМ время расчёта пропорционально N, поэтому соответствующее увеличение М приводит к тому, что время расчёта остаётся неизменным с ростом N.

1.1.3. Многомасштабные подходы

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

1.1.4. Стратегия развития МД modeling&simulation

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

В последнее время стало появляться значительное число экспериментальных работ, акцентирующих внимание на объёмном характере процесса модификации поверхности при лазерном облучении [7-9]. Однако расчёты, проводимые разными авторами, остаются или полностью одномерными (гидродинамическое моделирование), или квазиодномерными (атомистическое моделирование), когда только направление вглубь материала имеет микронный размер, а два

других направления имеют размеры в несколько нанометров и сшиваются через периодические граничные условия [10-17]. Поэтому целесообразным является переход к полномасштабному объёмному атомистическому моделированию процесса модификации, что является логическим развитием модели, но требует достаточно больших вычислительных ресурсов. Стоит отметить, что менее затратное гидродинамическое моделирование также может быть развито до объёмного моделирования, но для этого необходимо будет решить ряд принципиальных задач и, в частности, научиться корректно моделировать гидродинамику поверхностных волн. Однако даже при создании трёхмерной гидродинамической модели точность её будет значительно уступать атомистическому моделированию, так как при лазерном наноструктурировании большую роль играют такие процессы как плавление, разрушение и нуклеация, т.е. процессы, учёт которых в гидродинамическом подходе вызывает особую трудность. С другой стороны, в атомистическом моделировании все эти процессы описываются без привлечения каких-либо дополнительных приближений.

Остановимся более детально на модели, используемой в атомистическом моделировании для описания лазерного нанострук-туирования. При взаимодействии короткого лазерного импульса с веществом происходит сильный нагрев электронной подсистемы. Характерное время установления равновесия между электронами и ионами соизмеримо со временем самой модификации вещества и временами кинетических процессов, происходящих при этом (теп-лоперенос, фазовые переходы, возникновение ударных волн и т.д.) Таким образом, начальная стадия является двухтемпературной, когда температура электронов на порядок превышает температуру ионов. Для моделирования процесса лазерной модификации поверхности в цикле работ [10,13, 14,17] авторами развивается атомистическая модель двухтемпературного состояния. В данной модели используется приближение сплошной среды для электронной подсистемы и молекулярно-динамическое моделирование для ионной подсистемы. Таким образом, совместно решается система уравнений Ньютона для ионов и кинетическое уравнение теплопроводности для электронной подсистемы. Одной из особенностей развиваемой модели является учёт влияния электронного давления на динамику ионов.

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

На рис. 5 различные приближения атомистического моделирования процесса модификации поверхности схематично сопоставлены с числом вычислительных ячеек М (фактически с числом вычислительных ядер) и числом атомов Ж, требуемых для их реализации. Приближением, требуемым меньше всего вычислительных ресурсов, является квазиодномерный расчёт (1 на рис.5), когда вычислительная ячейка имеет размер около микрометра только в одном направлении (вглубь вещества), а два других направления имеют размеры в несколько нанометров и сшиваются через периодические граничные условия. Именно на таком типе задач отрабатывается и тестируется теоретическая основа модели. Стоит отметить, что необходимость микронного размера в направлении вглубь вещества обусловлена ударно-волновой природой лазерной абляции вещества и тем, что фазовые и структурные превращения происходят в поверхностном слое толщиной в несколько сотен нанометров.

Вторым по вычислительным затратам приближением является также квазиодномерный расчёт (2 на рис.5), но уже с размерами в несколько десятков нанометров по двум другим направлениям.

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

Третьим приближением является квазидвумерный расчёт (3 на рис.5), когда одно из направлений вдоль облучаемой поверхности имеет размер в несколько сотен нанометров. Такое приближение позволяет учитывать профиль энерговклада лазерного импульса. Стоит отметить, что в первых двух приближениях подразумевалось, что облучение происходит по всей поверхности с одинаковой интенсивностью, что конечно не соответствует действительности. Квазидвумерное приближение позволяет ввести распределение интенсивности по поверхности, но в двухмерной геометрии. И хотя реальный диаметр лазерного пятна в эксперименте составляет несколько микрометров, однако уже такое приближение является существенным достижением по сравнению с первыми двумя. Для непосредственного сравнения с экспериментом нужно использовать размерный анализ и теорию подобия, чтобы экстраполировать результаты моделирования до реальной ширины распределения интенсивности по поверхности.

Четвёртое приближение является также квазидвумерным (4 на рис.5), но уже с микронным размером в двух направлениях. Такая модель позволяет рассмотреть лазерное наноструктурирование за счёт плавления вещества. По результатам проведённых расчётов было установлено, что при небольших энерговкладах модификация поверхности может происходить не только за счёт лазерной абляции, но и за счёт плавления и расплескивания металла (прямая аналогия с капиллярными волнами на поверхности). Данный процесс хорошо демонстрирует, что при переходе к моделированию на больших масштабах существует возможность обнаружить новые коллективные явления атомов, недоступные для моделирования на меньших масштабах. Выполненные расчёты хорошо согласуются с недавними экспериментами (см. например [8,14]), где обнаружена модификация поверхности при низких энерговкладах в отсутствии абляции. Кроме того, моделирование в третьем и четвертом приближениях подтвердило высказанную ранее гипотезу, что размер зоны наноструктурирования определяется плавлением. В данной работе был сформулирован критерий, позволяющий рассчитать диаметр

зоны модификации во всем диапазоне энерговкладов. Граница зоны модификации поверхности, облучаемой лазерным импульсом, определяется той точкой пространства, где локальный энерговклад соответствует пороговому энерговкладу плавления. На рис. 6 показаны зависимости среднего радиуса поверхностной модификации от энергии импульса для золота и алюминия. Параметры импульса в расчётах выбирались из соответствия экспериментальному лазерному импульсу в работах [8,18] (соответственно, 1 и 2 на рис.6). Средний радиус определялся как среднее расстояние от центра лазерного пятна, где локальный энерговклад максимален, до точки, где локальный энерговклад соответствует пороговому энерговкладу плавления. Результаты расчётов на рис. 6 показаны линиями 3 и 4.

Пятое приближение подразумевает проведение полностью трёхмерного расчёта (5 на рис.5), на сопоставимых с экспериментом масштабах. Несмотря на то, что все остальные приближения могут дать много полезной информации о процессе лазерного нанострук-турирования, полностью описать динамику этого процесса возможно только в полномасштабном расчёте с микронными размерами во всех направлениях. Данный тип расчёта пока не был проведён, однако прогресс в вычислительных ресурсах, который наблюдается в последние годы, позволяет говорить о принципиальной возможности такого моделирования [18-20]. Проведение такого расчёта позволит напрямую сравнить эксперимент и атомистическое моделирование, опирающееся только на теорию из первых принципов. Такое сравнение важно не только само по себе, но и как возможность доказать или опровергнуть все предположения и гипотезы, которые были выдвинуты при моделировании на меньших масштабах.

1.1.5. Основания МД modeling & simulation

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

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

20 40 60 80

энергия лазерного импульса (нДж)

Рис. 6. Зависимости среднего радиуса поверхностной модификации от энергии лазерного импульса для золота и алюминия

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

Способы оценки значения N достаточно универсальны. Выбор его определяется масштабами пространственных и временных корреляций, характерных для поставленной задачи. В системе существует иерархия корреляций гс1 < гс2 < гс3 < ..., которой соответствует иерархия Nl < N2 < N3 < ..., где Ni = nr^i. Здесь n — концентрация частиц, rc — область расстояний, которая исследуется. Иными словами, можно сказать, что, выбирая то или иное значение N, мы тем самым обрываем ряд корреляций, которые можно будет исследовать в данном МД-расчёте. Выбирая N = nL3, мы также ограничиваем длины волн \ < L равновесных флуктуаций, т.е. фиксируем диапазон волновых векторов, для которого можно будет рассчитать дисперсию колебаний плотности: фононов в конденсированных средах, плазменных волн в неидеальной плазме, колебаний биомолекул и т.п. Подобным же образом выбор N ограничивает область исследуемых характеристик таких кооперативных явлений, как нуклеация, образование дислокаций и трещин и др.

Имеется иерархия времён корреляций: тс\ < тс2 < тс3 < .... Выбор Ь обрывает этот ряд, в частности, двумя неравенствами:

6< Ь2, а8тСг < Ь.

Здесь Б — коэффициент диффузии, а а8 — скорость звука.

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

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

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

1.2. Другие примеры масштабируемых задач

1.2.1. Радиационные повреждения в ядерных топливах

Физические модели, используемые в существующих топливных кодах, как правило, содержат широкий набор подгоночных параметров и эмпирических корреляций, что существенно снижает предсказательные возможности этих кодов, а также их применимость для описания новых видов ядерного топлива (в частности, нитридного и карбидного). В связи с этим, во всем мире предпринимаются попытки усовершенствования топливных кодов и создания механистических кодов нового поколения [21,22], где модели базируются на реальных физических механизмах (вместо упрощённых «корреляционных» моделей).

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

Используя метод согласования по силам («force-matching»), авторы разрабатывают новые многочастичные межатомные потенциалы для ураносодержащих материалов топлив для дальнейшего использования в МД-расчётах базовых микроскопических параметров моделей более высокого уровня [23-25]. Мировая практика создания потенциалов межчастичного взаимодействия показывает, что использования только данных экспериментов при разработке потенциалов оказывается недостаточно. Это связано в основном с тем, что набор экспериментальных данных, которые напрямую можно использовать при разработке потенциалов, зачастую сильно ограничен. В то же время, использование ab initio методов расчёта (таких как ТФП) сильно ограничено размерами систем из 1001000 атомов и очень короткими временами расчётов (для квантовой МД), чего зачастую недостаточно для расчёта поведения дефектов в решетке. Поэтому более подходящим на сегодняшний день вариантом является проведение расчётов в рамках классической МД, но с использованием потенциалов, описывающих свойства кристаллической решётки и энергетические характеристики дефектов в согласии с ТФП-расчётами (с привлечением экспериментальных данных, где они доступны).

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

На рис. 7 представлена перспектива расширения круга доступных процессов образования дефектов при увеличении числа ядер М в атомистических расчётах радиационных повреждений: 1 — МД и БЕТ расчёты энергий образования и миграции точечных дефектов; 2 — взаимодействие точечных дефектов друг с другом и кластерами, дислокационными петлями; 3 — генерация дефектов в столкновительных каскадах; 4—расчёты подвижности дислокационных петель и пор, их взаимодействий; 5—формирование дефектов при пролёте осколков деления в объёме и на поверхности материала (на данный момент такие расчёты ещё не проведены). Фрагменты моделирования и эксперимент взяты из работ авторов [21-27].

Определение механизмов и расчёты коэффициентов диффузии точечных дефектов являются наименее ресурсоёмкими среди прочих процессов радиационного повреждения, доступных для моделирования в рамках классической МД (уровень 1 на рис.7). Созданные потенциалы использованы для выяснения особенностей развития столкновительных каскадов в объёме и вблизи открытой поверхности материала (уровни 2-3 на рис. 7) и исследования положения линии плавления и объяснения расхождений в её экспериментальном определении на основе динамических и статических измерений.

Первостепенную важность для механистических моделей накопления радиационных повреждений и выхода продуктов деления представляют коэффициенты диффузии точечных дефектов различного типа и их небольших кластеров (в том числе вакансионных, например дефекта Шоттки в и02, и междоузельных дислокационных петель). Такие расчёты проведены для и02 [26], сплавов и-Мо [25]. Расчётные ячейки существенно большего размера (порядка 104-106 частиц) требуются для расчёта констант скоростей рекомбинации и кластеризации, эффективного радиуса взаимодействия точечных дефектов и мощности их стоков в поры, на поверхности [21,22]. Так в работе [21] проведён расчёт (для Мо) зависимостей от температуры ряда кинетических коэффициентов, необходимых для описания кинетики образования и роста дислокационных петель, полостей и пузырей Хе.

Расчёты генерации дефектов при облучении нейтронами, ионами, осколками деления весьма разнообразны (уровни 3-5 на рис. 7) из-за большого диапазона локальных энерговкладов. Доступными на сегодняшний день являются комбинированные модели [27], где пролёт высокоэнергетичного осколка деления многократно «разыгрывается» в рамках кинетического Монте-Карло, в результате получается распределение в пространстве и по энергии первичных выбитых атомов. А столкновительные каскады с энергией до ~ 100 кэВ, порождаемые этими первичными выбитыми атомами, рассчитываются на основе МД [23,27], в том числе с привлечением двухтемпературных моделей, учитывающих быструю ионизацию атомов и возбуждение электронной подсистемы. Результатом таких расчётов являются количество образовавшихся точечных дефектов, доля вакансионных и межузельных кластеров, их распределение по размерам и в пространстве (что может быть важно при учёте корреляций в последующих диффузионно-лимитируемых реакциях взаимодействия дефектов).

1.2.2. Нуклеация и кинетика фазовых переходов

Разрушение материалов при интенсивных импульсных воздействиях протекает через образование и распад метастабильных фазовых состояний. Распад перегретой или переохлаждённой жидкости, растянутого кристалла или жидкости происходит через образование зародышей новой фазы—нуклеацию. В большинстве моделей распада метастабильного состояния используется классическая теория нуклеации (КТН), разработанная в 1920-1950-е гг. в предположении

Рис. 7. Перспектива расширения круга доступных процессов образования дефектов при увеличении числа ядер М в атомистических расчётах радиационных поврежде-

слабого отличия свойств зародыша от свойств объёмной среды. В ряде работ указывается на то, что это предположение может не выполняться [28-30], вследствие чего КТН иногда даёт катастрофическое (5-10 порядков величины) несоответствие результатам экспериментов или прямого атомистического моделирования.

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

Рис. 8. Расширение круга доступных явлений при увеличении числа ядер для расчётов нуклеации полостей в растянутой жидкости. (изображение б) взято с сайта http://www.cse-lab.ethz.ch/)

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

Первой, наименее трудоёмкой задачей является наблюдение нуклеации единичной полости в заданном объёме (1 на рис. 8). Пространственный масштаб задачи определяется размером критического зародыша при заданной степени метастабильности, а временной — временем ожидания критического зародыша. Сложность задачи,

таким образом, растёт при движении от границы устойчивости фазы к кривой равновесия в силу роста размера критического зародыша и резкого увеличения времени жизни метастабильной фазы [28,31-33].

Второй этап—расчёт частоты нуклеации. Она является динамической характеристикой системы, и для её расчёта требуется усреднение времени ожидания зародыша по ансамблю для МД-запусков систем, находящихся в одном и том же термодинамическом состоянии (2 на рис.8). Эта задача допускает два способа распараллеливания [31-33]. Во-первых, каждый из расчётов в ансамбле может быть распараллелен при помощи декомпозиции по пространству. Во-вторых, для расчётов средних величин по ансамблю можно применить «декомпозицию по ансамблю», то есть параллельный расчёт нескольких независимых МД-траекторий на различных узлах многоядерной системы. В силу независимости МД-траекто-рий в ансамбле при втором виде распараллеливания достигается идеальная масштабируемость из-за отсутствия межпроцессорного взаимодействия.

Третьей по вычислительной сложности задачей является исследование роста зародыша новой фазы (3 на рис.8). Динамика роста полости в жидкости в кубической ячейке с ребром Ь может быть исследована лишь на временах, меньших ¿таж = Ь/с3, где с3 — скорость звука в жидкости, т.к. затем на рост полости существенно влияют выбранные граничные условия.

Четвёртый масштабный уровень возникает при дальнейшем увеличении размера моделируемой системы, когда достигается режим одновременного зарождения и роста множества (102-103) полостей в объёме расчётной ячейки при импульсном растяжении (4 на рис.8). Такие расчёты позволяют рассчитать предельно достижимое напряжение в микрообъёме жидкости при его однородном нагружении.

Пятый масштабный уровень может быть достигнут на мощнейших на данный момент суперкомпьютерах. Прямое МД-моделиро-вание объёма 1 мкм3 и более уже позволяет исследовать с атомным разрешением процессы в системе с неоднородно распределёнными термодинамическими характеристиками: эффекты прохождения ударной волны, её отражения от свободной поверхности и формирования откольной пластины (5 на рис.8). Отметим, что подобные

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

Для перехода к кинетике фазового перехода на мезо- и макромасштабах по данным МД-расчётов строятся кинетические модели фазового перехода. На основе рассчитанных частот нуклеации и скорости роста зародышей строятся аппроксимационные зависимости этих скоростей от давления и температуры. На основе этих зависимостей с помощью кинетической модели типа «нуклеация и рост» рассчитывается суммарный объём полостей и их распределение по размерам в произвольный момент времени при заданной истории нагружения [28,34]. На врезке а) к рис. 8 показано хорошее согласие распределения полостей по размерам, рассчитанного по многомасштабной модели, с полученным в прямом МД-моделиро-вании при той же скорости растяжения. Такой подход позволяет рассчитать откольную прочность (предел достигаемых в жидкости напряжений), а также время жизни жидкости до разрыва.

Введение кинетической модели фазового перехода в расчёт методом механики сплошной среды позволяет более точно рассчитывать результат импульсных механических и тепловых воздействий на вещество, в частности, рассчитывать скорость свободной поверхности и профили давления в образцах при ударно-волновом воздействии, предсказать микроструктуру откольной пластины, воздействие схло-пывания полостей (кавитации) на находящиеся в жидкости тела. В статье [35] приведены первые результаты подобного расчёта для откола в жидком алюминии. В [36] приводятся результаты рекордного гидродинамического расчёта схлопывания 15000 кавита-ционных полостей в воде (см. б) на рис.8). Отметим, что параметры нуклеации и роста зародышей для модели фазового перехода, а также уравнение состояния и коэффициенты переноса для модели сплошной среды могут быть рассчитаны методами классической или квантовой МД [28, 34,37]. Таким образом, возможно построение моделей поведения вещества при импульсном воздействии без необходимости введения свободных параметров.

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

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

1.2.3. Супрамолекулярные системы

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

Квантовый уровень используется в первую очередь для расчёта энергий взаимодействия молекул и создания потенциалов взаимодействия. Полученные результаты могут использоваться напрямую либо в термодинамических моделях, либо в методах молекулярной динамики и Монте-Карло. Термодинамические модели основаны на уравнениях статистической физики, записанных с теми или иными упрощениями, и позволяют получать фазовые диаграммы различных газовых гидратов и их степени заполнения. Такой способ даёт предсказания, достаточные для инженерной точности, однако не подходит для изучения свойств на микроуровне. Для этих целей лучше подходят методы молекулярной динамики и Монте-Карло. Они используются для изучения устойчивости структур, нуклеации и гидратообразования, процессов ингибирования, поверхностных, спектральных (получение колебательных спектров) и транспортных (диффузия, теплопроводность) свойств. Точность таких методов зависит от выбранного потенциала взаимодействия и требует дополнительной верификации. В работах авторов [38,39] приведены примеры использования этих методов для моделирования фазовой диаграммы и процессов нуклеации в гидратах метана, а также определения стабильности новых структур гидратов водорода.

число вычислительных ядер 101 102 103 104 105 10® I_I_I_I_I_I

число атомов

I_I_I_I_I_I_I_I

ю2 ю3 ю4 ю5 ю6 ю7 ю8 ю9

Рис. 9. Расширение круга доступных явлений при моделировании газовых гидратов

Для моделирования газовых гидратов используются преимущественно дальнодействующие потенциалы взаимодействия, что сильно снижает пространственные и временные размеры моделируемой системы. Обычный размер расчётной ячейки не превышает десятков тысяч атомов, что может приводить к появлению размерных эффектов и затрудняет прямое сравнение с экспериментом. Дальнейшее развитие связано с распространением суперкомпьютеров типа Blue Gene с тороидальной архитектурой, что позволит уже в ближайшей перспективе вывести расчёты на уровень миллионов атомов. Возможность эффективного распараллеливания МД-расчётов газовых гидратов на Blue Gene продемонстрирована в статье [40]. Дальнейший переход на эксафлопсный уровень суперкомпьютеров позволит моделировать микрокристаллы с ярко выраженными размерными эффектами. Возможность такого перехода и перспектива дальнейшего развития схематично представлена на рис.9: от термодинамической устойчивости (1) и трёхфазного

моделирования плавления (2) к исследованию кинетического инги-бирования роста (3) и образованию гидрата на поверхности льда (гетерогенная нуклеация, 4).

1.2.4. Полимеры

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

Свойства полимерных нанокомпозитов существенным и не вполне очевидным образом зависят от целого ряда параметров: типа, концентрации, ориентации, упругих свойств наполнителя, пластических и термических свойств самой полимерной матрицы. Описание подобной многопараметрической задачи на микроуровне с использованием сугубо экспериментальных методов крайне затруднительно и пока во многом носит эмпирический характер. В то же время предсказательные возможности методов компьютерного моделирования достаточно давно достигли уровня, позволяющего как изучать природу химических процессов, происходящих на атомарном уровне, так и производить качественную и количественную оценку макропараметров нанокомпозитных материалов. Если в первом случае речь идёт о квантовомеханических и классических полноатомных молекулярно-динамических расчётах, то вторая ситуация требует привлечения методов, позволяющих выйти за пространственные (и временные) масштабы, характерные для типичных МД-расчётов. На помощь в этом случае приходят методы огрублённой (т.н. coarse-grained) МД.

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

finite element

RnfiMfl

микросекунды -

наносекунды

пикосекунды

фемтосекунды —

1 нм

1-10 нм 10-1000 нм 1-100 им 1-10 мм

Рис. 10. Принципиальная схема многомасштабного моделирования полимерных композитов

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

Связь масштабов помимо сохранения характера межмолекулярных взаимодействий обеспечивается передачей параметров структуры материала между уровнями. Ценность подобного последовательного с физической точки зрения подхода определяется тем, что основанные на нем многоуровневые модели могут иметь не только описательную, но и прогностическую силу, что чрезвычайно важно для решения задачи о создании материала с заданными свойствами и для разработки его состава in silico (рис.10). Представленную на иллюстрации иерархию элементарных процессов необходимо прослеживать, начиная с атомистического уровня классических и квантовых моделей межатомного и межмолекулярного взаимодействия, т.к. только в этом случае создаваемые многомасштабные модели будут способны учитывать нюансы изменения свойств материала при изменении деталей его атомарной структуры. В качестве иллюстрации ускорения вычислений при переходе к огрублённой системе и проверки её масштабируемости был проведён ряд расчётов с моделью полиэтиленовой матрицы, насчитывавшей 30000

Тпс (секунды)

1 -1-1-Г-1-1-1—|—Г^-1

100 1000 Мощность задействованного вычислительного ПОЛЯ (Гфлопс)

Рис. 11. Зависимость вычислительного времени, требуемого на расчёт траектории длительностью в одну пикосекунду (Тпс), от величины задействованных вычислительных ресурсов

атомов. На рис.11 представлена зависимость реального вычислительного времени, требуемого на расчёт траектории длительностью в одну пикосекунду (Тпс), от величины задействованных вычислительных ресурсов. Расчёт производился в пакете LAMMPS с использованием потенциала AIREBO [41] для полноатомной модели и DREIDING [42] для огрублённой. Силовые константы последнего предварительно вычислялись с помощью методики согласования по силе (force-matching method) на основе поведения полноатомной модели. В качестве супер-атомов рассматривались CH2 группы.

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

огрублённой модели обеспечил ускорение расчёта более чем в 30 раз.

Концентрации наполнителя в нанокомпозитах невелики, их массовые доли, как правило, не превышают 4-5%. Помимо того, что сами нановключения в контексте метода молекулярной динамики являются достаточно крупными объектами, содержащими от нескольких десятков тысяч атомов, для расчёта макропараметров материала требуется задействовать в модели существенное их количество. Таким образом, задача анализа свойств полимерных нанокомпозитов предполагает проведение расчётов со сверхкрупными системами, отражающих при этом тонкости взаимодействия полимера с нановключением на атомарном уровне. Многомасштабный подход в подобных задачах потенциально позволяет выйти к системам с размерами порядка десятков микрометров.

1.2.5. Механика сплошных сред

Картинки, аналогичные вышеприведённым, могут быть построены и для задач, решаемых численными методами МСС. По нашему предложению был построен пример для задач сейсмикии [43].

1.3. Кулоновские системы

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

1.3.1. Пылевая плазма

Рассмотрим данный подход на примере моделирования пылевой плазмы [44]. Пылевая плазма является ионизованным газом, содержащим макрочастицы конденсированного вещества (пыль) микронных размеров. Такие частицы либо самопроизвольно образуются в плазме в результате различных процессов, либо вводятся в плазму извне. Пылевая плазма встречается в хвостах комет, кольцах планет-гигантов, верхних слоях атмосфер планет, и может образовываться в результате технологических процессов (при сгорании топлив, травлении и напылении, в производстве наночастиц, при работе токамака и т.д.). Наличие макрочастиц в плазме может

существенно влиять на её свойства, а также на термодинамику системы и процессы переноса. Макрочастицы, изначально не обладающие ни электрическим зарядом, ни дипольным (квадру-польным и т.д.) моментом, в плазме могут заряжаться потоками электронов и ионов, а также путём фото-, термо- или вторичной эмиссии электронов и приобретать значительный отрицательный или положительный электрический заряд. Большинство исследуемых явлений в пылевой плазме имеет синергетическую природу из-за самосогласованности системы пылевых частиц и плазмы, а также из-за сильного взаимодействия пылевых частиц друг с другом. Теоретическое исследование плазменно-пылевой системы и её теплофизических свойств основывается на методах статистической физики, классической механики и МД-моделирования [45-49]. Учитывается анизотропия температуры системы пылевых частиц.

В рамках рассматриваемых подходов было осуществлено предварительное исследование условий применимости термодинамических функций к пылевой плазме [49]. Создана методика по исследованию свойств плазменно-пылевой системы в газовом разряде с помощью метода МД, в частности, написаны уравнения движения пылевых частиц. Исследовано явление аномального разогрева колебаний пылевых частиц в плазме и предложен механизм разогрева [45-48]. Обнаружена возможность резонанса колебаний пылевых частиц, откуда сделаны оценки частоты, амплитуды и кинетической энергии колебаний пылинок. Написан программный пакет, моделирующий движение системы пылевых частиц в газовом разряде при комнатной температуре с учётом флуктуаций заряда и особенностей приэлектродного слоя. Наблюдаемая в эксперименте анизотропия кинетической температуры пылевых частиц была обнаружена и в моделировании.

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

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

На рис. 12 представлены приближения, используемые при моделировании плазменно-пылевой системы, и качественно сопоставлены с числом необходимых вычислительных ядер и числом пылевых частиц в расчёте. Моделировать одновременно и плазму, и большое число пылевых частиц на данный момент не представляется возможным, поэтому используется подход, при котором пылевая плазма моделируется поэтапно. Сначала исследуется влияние плазмы на пылевую частицу (1 на рис. 12), процессы зарядки пылевой частицы и особенности флуктуаций заряда пылевой частицы. В частности, для исследования динамики пылевых частиц важное значение имеет учёт в модели зависимости флуктуаций заряда от флуктуаций потоков электронов и ионов, от движения самой пылевой частицы и соседних пылевых частиц, от влияния приэлек-тродного слоя газового разряда, от флуктуаций движения пылевых частиц из-за окружающего газа. На втором этапе (2 на рис 12) исследуются перенос энергии между степенями свободы, зависимость характеристик системы от её параметров, механизмы разогрева и процесс релаксации системы. Вычислительная сложность этого этапа обусловлена резким отличием временных масштабов характерных процессов в системе, что требует больше миллиарда шагов для вывода системы на квазиравновесие. Необходимо проводить варьирование 15 параметров, чтобы заполнить сетку в 15-мерном пространстве, где каждый узел сетки соответствует какому-то набору параметров и, как минимум, одному численному эксперименту. На этом этапе определяются термодинамические коэффициенты, области их применения, связи между параметрами системы. Уже на данном этапе проводится частичное сравнение с экспериментом.

Третий этап (3 на рис 12, [44]) посвящён изучению свойств облака пылевых частиц в плазме, в частности, пылезвуковых волн, особенностей поведения больших структур пылевых частиц разных размеров. В данном приближении используются коэффициенты, полученные на втором этапе. Первые два этапа дают значительную

• •

10° ю1 ю2 ю3 ю4 '-'-'-'-'-' N

10° ю1 ю2 ю3 ю4 ю5

Рис. 12. Этапы моделирования пылевой плазмы в различных приближениях в зависимости от числа вычислительных ядер М и рассматриваемого числа пылевых частиц N. У каждой картинки представлены характерные размеры моделируемой области

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

1.3.2. Ионная плазма и пробой газов и жидкостей

При импульсных разрядах в плотных электроотрицательных средах в результате прилипания электронов за времена порядка наносекунд образуется ионная плазма. Характерное время жизни

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

В работах [50-53] рассмотрено два фактора подавления рекомбинации: «узкое горло» и сольватация ионов. Развиты подходы, сочетающие молекулярное моделирование и аналитическое рассмотрение.

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

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

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

Рис. 13. Масштабируемость моделирования ионных систем. 1 — модель ионной системы с эффективным учётом нейтральной компоненты; 2 — модель ионной системы с детальным описанием нейтральной компоненты в рамках единой модели; 3 — методы макроскопической кинетики

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

Коэффициенты скоростей элементарных процессов, полученные в результате моделирования 1 или 2, используются в методах макроскопической кинетики 3 для описания вещества в изучаемой системе.

м

10° 102 103 104 105 106

Рис.14. Масштабируемость задач сольватации. 1 — вычисление потенциалов взаимодействия ионов и молекул методами квантовой химии; 2 — классическая модель кластера; 3—макроскопические модели

Сольватация ионов. Рекомбинация ионов в плазме фтора замедлена сольватацией ионов и проходит через образование промежуточной пары кластерных ионов. В этом случае необходим детальный учёт формирования ион-молекулярных кластеров. Масштабируемость задач схематично представлена на рис. 14.

Этап 1 включает расчёт энергий взаимодействия нескольких рассматриваемых частиц методами квантовой химии. На основании полученных в таком расчёте потенциалов межчастичного взаимодействия на этапе 2 осуществляется вычисление структур, энергий и взаимодействия кластеров, формируемых из рассматриваемых ионов и молекул в рамках классических молекулярно-динамических подходов. Результатом этих расчётов являются данные, которые позволяют вычислить поправки к константам скоростей элементарных процессов с участием кластерных ионов. Полученные кинетические коэффициенты используются на этапе 3 для макроскопического кинетического и гидродинамического моделирования системы.

1.3.3. Электрохимия и электрический двойной слой

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

Двойной электрический слой в суперконденсаторах образуется на границе проводника и жидкого или твёрдого вещества, обладающего ионной проводимостью. Проводимость проводника может иметь электронную, дырочную или электронно-дырочную природу. В соответствии с типом используемого электролита устройства подразделяются на несколько групп. К первой группе относятся конденсаторы с твёрдым электролитом, в качестве которого чаще всего находят применение двойные соли серебра, обладающие высокой ионной проводимостью при комнатной температуре. Такие конденсаторы, как правило, имеют сравнительно небольшие рабочее напряжение (порядка 0.7 В) и энергоёмкость. Это компактные твердотельные устройства, находящие применение в микроэлектронике. Ко второй группе относятся конденсаторы с водным электролитом, представляющим чаще всего водный раствор щелочей или серной кислоты. Рабочее напряжение в таких системах может составлять величину до 1.2 В, а энергоёмкость до 100 Дж/г. Такие конденсаторы обладают низким внутренним сопротивлением и высокой удельной мощностью, что объясняется высокой подвижностью протона в водных растворах. Наибольшими значениями удельной энергоёмкости обладают конденсаторы третьей группы, в которых в качестве электролита используются апротонные растворители, как правило, ионные жидкости. Конденсаторы с неводными электролитами обычно уступают конденсаторам с кислотными и щелочными электролитами по величине удельной мощности. Независимо от типа используемого электролита, в качестве материала электродов в суперконденсаторах используются активированный уголь, углеродные нанотрубки, графен.

В работе было осуществлено предварительное исследование двойного электрического слоя на границе плоского слоя графита и водного раствора щелочного электролита [54]. Получены следующие результаты:

(1) влияние двойного электрического слоя в электронно-дырочной плазме графита на ёмкость системы преобладает по сравнению с влиянием двойного слоя в электролите на плоском слое графита;

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

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

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

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

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

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

10° ю2 103 104 105 10е ю7

Рис. 15. Расширение круга явлений и задач с увеличением числа ядер при моделировании двойного слоя. 1 — зонная структура и экранирование в материале электрода; 2 — квантовая модель нанопоры; 3 — потенциал взаимодействия молекул с поверхностью; 4 — классическая модель микропоры; 5 — модель переноса ионов в процессе заряда и разряда в угле с реалистичной топологией пор

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

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

Рис. 16. Сопоставление масштабируемости на системах толстого дерева (слева) и тора Б1ибСбПб/Р, ЛМЬ (справа) для леннард-джонсовской системы. В обоих случаях использован пакет ЬЛММРЯ

разрядки в углях с реалистичной топологией пор, что является исключительно ресурсоёмкой задачей — этап 5.

2. Какие нужны суперкомпьютеры эксафлопсного класса?

Рассмотрим три класса задач.

2.1. Классическое атомистическое моделирование

Результаты, представленные на рис. 16, показывают, что, например, для 108 частиц снижение эффективности распараллеливания до 80% для толстого дерева происходит при 700 процессорах, когда ускорение составляет всего « 600 раз, в то время как для тора—при «120000 процессорах, когда ускорение составляет уже « 105 раз. Преимущества тороидальной системы весьма значительны. В силу этого, соотношения между осями М и N на рисунках, приведённых в разделе 1.2, будут различаться для систем толстого дерева и тора.

2.2. Механика сплошных сред

Ограничение масштабируемости алгоритмов как молекулярной динамики, так и механики сплошных сред (МСС), определяется межпроцессорным обменом данными, что схематично проиллюстрировано на рис. 2 и 18. При декомпозиции по пространству система

128 512 2048 8192 32768 131072

Рис. 17. Распределение потерь времени для тороидальной системы

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

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

В [56] представлен новый метод решения задач вычислительной гидродинамики на суперкомпьютерах петафлопсной производительности. Продемонстрирована масштабируемость, близкая к идеальной, вплоть до 100 000 процессоров на BlueGene/P(ANL). В [43] продемонстрирована масштабируемость задач сейсмики, близкую к идеальной, вплоть до 10 000 процессоров.

Рис. 18. Схематичное сопоставление масштабируемости атомистического моделирования (слева, см. рис. 2) и в задачах механики сплошных сред (справа)

2.3. Квантовое (ab initio) атомистическое моделирование

2.3.1. Параллелизм и масштабируемость

Хотя наличие более глубокого, квантового уровня подразумевалось и в предыдущих разделах, его масштабируемость дана весьма условно, поскольку по сравнению с одномерной масштабируемостью (рис.3) появляется ещё одна ось — число базисных функций, по которым раскладывается квантовое решение (рис. 19). В работе [57] представлен обзор распараллеливания алгоритмов ab initio МД на основе теории функционала электронной плотности в базисе плоских волн. Проанализированы требования к балансу вычислительной мощности узлов и коммутационной сети суперЭВМ с точки зрения достижения максимальной эффективности для примеров особенно требовательных в вычислительном отношении задач физики разогретого плотного вещества. Описана альтернативная стратегия параллелизации в вейвлетном базисе и выигрыш в производительности при использовании гибридных вычислительных систем в этом случае.

2.3.2. Роль архитектуры

Как и МД-задачи, ab initio расчёты являются чувствительными к архитектуре вычислительных кластеров. Вклад в снижение масштабируемости вносит обмен данными между узлами, поскольку

10

Рис. 19. Иллюстрация многопараметрической масштабируемости задач квантового моделирования

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

На рис. 20 представлено сравнение эффективности масштабирования ТФП-расчёта для двух суперкомпьютеров с различной топологией коммутационных сетей. В качестве тестового примера производительности был выбран пакет CP2K. Расчёты проводились для системы 2048 молекул воды. Видно, что на системе с топологией толстое дерево данный тип задач в пределах больших размеров вычислительного поля масштабируется хуже: существенное отклонение от прямой идеального масштабирования наблюдается уже при использовании 512 ядер, в то время как на IBM BlueGene/P [58] ускорение растёт линейно вплоть до подключения нескольких тысяч ядер.

1Д

ч—I—ивм вар

О 10000 20000 зоооо

1 1 I ' I 1 I ' Число ядер К-100

0 200 4СС 600 800 1000

Число ядер Вй/Р

1-1-1-1-1-1-1-1-1-1-1

О 2000 4 С 00 6000 8000 10000

Рис. 20. Оценка параллельной эффективности для модели 2048 молекул воды для различных топологий суперкомпьютеров (расчёт в рамках ТФП с использованием пакета CP2K). Приведены данные для топологии толстое дерево (суперкомпьютер К-100) и 3-х мерный тор

2.4. Локальность обменов данными

Нелокальные обмены, однако, являются распространённой чертой вычислительных алгоритмов. Широкий класс подобных алгоритмов в области молекулярного моделирования соответствует многомерному быстрому преобразованию Фурье. Адаптация подобных алгоритмов на сети с тороидальной топологией является уже весьма разработанной частью необходимого инструментария параллельных вычислительных алгоритмов [59].

Тороидальные топологии представляют одно из магистральных направлений развития интерконнекта. На топологии 5-ти мерного тора основаны суперкомпьютеры IBM BlueGene/Q, неполному 6-ти мерному тору соответствует интерконнект Fujitsu Tofu. Таким образом, развитие идёт в направлении роста локальной связности коммуникаций, хотя «тороидальность» сама по себе и не является обязательной. Так на смену интерконнекту с топологией 3-х мерного тора Cray Gemini пришла новая коммуникационная сеть Cray Aries с гибридной топологией «стрекоза» (dragonfly).

Очевидно, что развиваемые вычислительные алгоритмы «экса-флопсной эры» должны будут учитывать особенности интерконнекта. И наоборот: выбор коммуникационной сети должен определяться будущими задачами. При этом проблема наилучшего отображения коммуникационного шаблона задачи на топологию коммуникационной сети выходит на передний план [60] вместе с требованиями повышения локальности коммуникаций используемых алгоритмов для снижения энергозатрат [61].

2.5. Графические ускорители

Среди современных высокопроизводительных систем все чаще используются гибридные системы, включающие графические ускорители (ГУ) или иные специализированные устройства для выполнения векторных операций в рамках модели SIMD. Этому во многом способствовало появление сред программирования, таких как CUDA и OpenCL, предназначенных для создания и выполнения на ГУ программ, не связанных с обработкой изображений. Гибридные системы на основе ускорителей зачастую демонстрируют быстродействие, в десятки или сотни раз превышающее быстродействие систем на традиционных процессорах (ЦПУ—центральное процессорное устройство). Однако создание программ, эффективно использующих вычислительные возможности ускорителей, это достаточно трудоёмкий процесс, требующий учёта специфики аппаратной архитектуры, без которого производительность на конкретной задаче может оказаться существенно меньше теоретического предела для данного устройства.

Отличительными особенностями ГУ являются большое количество простых вычислительных ядер, аппаратная поддержка мно-гопоточности (multi-threading), высокая пропускная способностью внутренней памяти, оптимизированной под последовательный доступ, при сравнительно медленной шине данных, связывающей ГУ и основную память вычислительного узла. Это определяет требования к алгоритмам и набор задач, которые могут быть эффективно решены с применением ГУ. В первую очередь алгоритмы должны обладать высокой степенью параллелизма по данным, минимальным количествам обращений к памяти по произвольным адресам, локальностью данных, преобладанием вычислений над операциями ввода-вывода.

Алгоритмы МД во многом удовлетворяют перечисленным выше требованиям, поэтому за последние несколько лет появились программы МД-моделирования, позволяющих проводить расчёты с использованием ГУ. Это программы LAMMPS [62], HOOMD-blue [63], NAMD и VMD [64], GROMACS [65], HALMD [66], и др. Большинство из перечисленных программ является усовершенствованием существовавших ранее кодов, в то время как пакет HOOMD-blue изначально создавался для вычислений на ГУ, в результате чего его производительность превосходила, например, ранние реализации ГУ-алгоритмов в LAMMPS [67]. Как было указано ранее, алгоритмы МД отличаются высокой эффективностью распараллеливания на традиционных (ЦПУ) кластерах. Вместе с тем, следовало бы ожидать, что эффективность распараллеливания на гибридных кластерах, содержащих ГУ, окажется значительно меньше за счёт того, что при увеличении скорости вычислений на отдельном узле кластера относительный вклад потерь на передачу данных по внутренней сети (interconnect) окажется выше.

На рис. 21 показана скорость работы программы LAMMPS, т.е. количество МД-шагов в секунду, в зависимости от числа вычислительных элементов: процессорных ядер для ЦПУ и ускорителей для гибридной системы. В тестах моделировалась леннард-джон-соновская жидкость с температурой Т = 1.0, плотностью р = 0.19, отсечкой потенциала rcut = 3 (в единицах LJ). Расчёты проводились на гибридном кластере К-100 ИПМ РАН, узел которого включает 12 процессорных ядер Intel Xeon X5670 и 3 ускорителя NVidia Tesla C2050. Вопреки негативным ожиданиям, приведённые результаты показывают достаточно высокую эффективность распараллеливания: 62% для ГУ при 80% для ЦПУ. Это позволяет говорить

Рис. 21. Эффективность распараллеливания: скорость работы программы (числа МД-шагов в секунду) на ЦПУ (нижний ряд точек) и ГУ (верхний ряд точек) в зависимости от числа процессорных ядер для нижнего ряда или числа ускорителей для верхнего. Пунктирные линии — идеальное распараллеливание, сплошные — степенная аппроксимация; (а) 1 млн. частиц, (б) 4 млн. частиц в ячейке

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

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

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

Рис. 22. Характерные размеры МД-ячейки для моделирования систем с кулоновским взаимодействием и области применения традиционных процессоров (CPU) и графических ускорителей (GPU)

по сравнению с одним ядром ЦПУ может достигать 100 раз и позволяет на два порядка увеличить размер исследуемой системы [68]. При этом системы до 105 частиц могуть рассчитываться на одном узле кластера, а коммуникации между узлами применяться для лишь статистического усреднения. Дополнительное увеличение размеров моделируемой системы возможно за счёт использования приближенных методов, например «Particle-particle particle-mesh» [69].

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

3. Заключение

3.1. Зачем нужны суперкомпьютеры эксафлопсного класса

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

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

3.2. Какие нужны суперкомпьютеры эксафлопсного класса

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

Показано, что использование суперкомпьютеров с тороидальной топологией оказывается предпочтительным для всех трёх классов задач. При этом речь может идти о выигрыше в производительности в десять и более раз, а в некоторых случаях задача, решаемая на тороидальной системе, оказывается не решаемой на системе толстое дерево. Заметим, что суперкомпьютер «Млечный Путь» (№1 в Тор500) может провалиться вниз для конкретной научной задачи ввиду своей сомнительной архитектуры.

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

3.3. Взаимодействие разработчиков и пользователей суперкомпьютеров

В новой редакции Top500 [71] от 19 ноября 2013 в верхней десятке стало на 1 тороидальную систему больше. Было 6, стало 7. Это связано с тем, что суперкомпьютеры в США создаются под определённые классы задач с учётом мнения пользователей. Два примера из практики IBM разных десятилетий.

При разработке архитектуры суперкомпьютеров IBM BlueGene/L задачи молекулярного моделирования рассматривались в качестве главного приоритета [2].

А ещё в году 1970-ом один из авторов статьи (Г. Норман) участвовал в беседе с приехавшим в Москву Энрико Клементи, видным сотрудником IBM, разработчиком квантовохимической программы IBMOL. Группа Клементи публиковала тогда десятки статей в лучших журналах. На вопрос, на каких экземплярах новых машин IBM он считает, под номером один или два, Энрико небрежно ответил -- под номером Ноль.

Благодарности. Работа частично поддержана грантом РФФИ 13-01-12070-офи_м, а отдельные авторы—13-08-01022-а (Ланкин А. В., Норман Г. Э.), 13-08-01428-а (Куксин А.Ю., Смирнов Г. С., Стегайлов В. В.), 12-08-00666-а (Орехов Н.Д., Стариков С. А.), 12-08-33140-мол-а_вед (Куксин А. Ю., Писарев В. В., Стариков С. А.), 12-02-33170 мол-а_вед (Ланкин А. В., Морозов И. В.), 12-02-31783-мол-а (Морозов И. В.), 12-08-31461-мол_а (Тимофеев

A. В.), 14-08-31550-мол_а (Смирнов Г. С., Орехов Н.Д.), 14-08-31587-мол_а (Писарев В. В.), стипендией Президента РФ для молодых учёных (Писарев В. В.). Расчеты выполнены в МСЦ РАН, на кластерах МФТИ-60 кафедры информатики МФТИ, К-100 ИПМ им. М. В. Келдыша РАН, «Ломоносов» МГУ им. М.

B. Ломоносова. Использовались пакеты LAMMPS, VASP, CPMD, СР2К и собственные разработки, о которых сказано в тексте.

Описок литературы

[1] INCITE program, URL http://www.doeleadershipcomputing.org/. f192

[2] IBM Journal of Research and Development, 2005, URL http://www.research. ibm.com/journal/rd49-23.html. f192, 236

[3] В. В. Стегайлов, Г. Э. Норман. Проблемы развития суперкомпьютерной отрасли в России: взгляд пользователя высокопроизводительных систем // Программные системы: теория и приложения, 2014. Т. 5, № 1, с. 111— 152, URL http://psta.psiras.ru/psta2014_1_111-152.pdf. f!94

[4] Г. Э. Норман, В. В. Стегайлов. Стохастическая теория метода классической молекулярной динамики // Математическое моделирование, 2012. Т. 24, № 6, с. 3-44. f194

[5] А. В. Янилкин, П. А. Жиляев, А. Ю. Куксин, Г. Э. Норман, В. В. Писарев, В. В. Стегайлов. Применение суперкомпьютеров для молекулярно-динамического моделирования процессов в конденсированных средах. // Вычислительные методы и программирование, 2010. Т. 11, с. 111-116. ^196

[6] Z.A. Insepov, E. M. Karatajev, G. E. Norman. The kinetics of condensation behind the shock front // Zeitschrift fur Phys. D Atoms, Mol. Clust., 1991. Vol. 20, no. 1-4, p. 449-451. f197

[7] J. M. Savolainen, M. S. Christensen, P. Balling. Material swelling as the first step in the ablation of metals by ultrashort laser pulses // Phys. Rev. B, 2011. Vol. 84, p. 193410. f197

[8] M. Ishino, A. Y. Faenov, M. Tanaka, S. Tamotsu, N. Hasegawa, M. Nishikino, T. A. Pikuz, T. Kaihori, T. Kawachi. Observations of surface modifications induced by the multiple pulse irradiation using a soft picosecond x-ray laser beam // Appl. Phys. A, 2012. Vol. 110, p. 179-188. f197, 200, 201

[9] С. И. Ашитков, Н. А. Иногамов, В. В. Жаховс кий, Ю. Н. Эмиров, М. Б. Агранат, И. И. Олейник, С. И. Анисимов, В. Е. Фортов. Образование на-нополостей в поверхностном слое алюминиевой мишени при воздействии фемтосекундных лазерных импульсов // Письма в ЖЭТФ, 2012. Т. 95, № 4, с. 195-197. f197

[10] С. В. Стариков, В. В. Стегайлов, Г. Э. Норман, В. Е. Фортов, М. Ишино, М. Танака, Н. Хасегава, М. Нишикино, Т. Охба, Т. Каихори, Е. Очи, Т. Имазоно, Т. Кавачи, С. Тамотсу, Т. А. Пикуз, И. Ю. Скобелев, А. Я. Фаенов. Лазерная абляция .золота: эксперимент и атомистическое моделирование // Письма в ЖЭТФ, 2011. Т. 93, № 11, с. 719-725. f198

[11] M. E. Povarnitsyn, T.E. Itina, P. R. Levashov, K. V. Khishchenko. Simulation of ultrashort double-pulse laser ablation // Appl. Surf. Sci., 2011. Vol. 257, p. 5168-5171. f198

[12] A. V. Mazhukin, V. I. Mazhukin, M. M. Demin. Modeling of femtosecond ablation of aluminum film with single laser pulses // Appl. Surf. Sci., 2011. Vol. 257, p. 5443-5446. f198

[13] Г. Э. Норман, С. В. Стариков, В. В. Стегайлов. Атомистическое моделирование лазерной абляции золота: эффект релаксации давления // ЖЭТФ, 2012. Т. 141, №5, c. 910-918. f198

[14] G. Norman, S. Starikov, V. Stegailov, V. Fortov, I. Skobelev, T. Pikuz, A. Faenov, S. Tamotsu, Y. Kato, M. Ishino, M. Tanaka, N. Hasegawa, M. Nishikino, T. Ohba, T. Kaihori, Y. Ochi, T. Imazono, Y. Fukuda, M. Kando, T. Kawachi. Nanomodification of gold surface by picosecond soft x-ray laser pulse // J. Appl. Phys., 2012. Vol. 112, no. 1, p. 013104. f198, 200

[15] M. E. Povarnitsyn, T. E. Itina, P. R. Levashov, K. V. Khishchenko. Mechanisms of nanoparticle formation by ultra-short laser ablation of metals in liquid, environment // Phys. Chem. Chem. Phys., 2013. Vol. 15, p. 3108-3114. f198

[16] N. A. Inogamov, V. V. Zhakhovsky, Yu. V. Petrov, V. A. Khokhlov, S.I. Ashitkov, K.V. Khishchenko, K.P. Migdal, D. K. Ilnitsky, Yu. N. Emirov,

P. S. Komarov, V.V. Shepelev, C.W. Miller, I.I. Oleynik, M. B. Agranat, A. V. Andriyash, S. I. Anisimov, V. E. Fortov. Electron-ion relaxation, phase transitions, and surface nano-structuring produced by ultrashort laser pulses in metals // Contrib. to Plasma Phys., 2013. Vol. 53, no. 10, p. 796-810. f198

[17] G.E. Norman, S.V. Starikov, V.V. Stegailov, I.M. Saitov, P. A. Zhilyaev. Atomistic Modeling of Warm Dense Matter in the Two-Temperature State // Contrib. to Plasma Phys., 2013. Vol. 53, no. 2, p. 129-139. f198

[18] S.V. Starikov, A.Ya. Faenov, T. A. Pikuz, I.Yu. Skobelev, V. E. Fortov, S. Tamotsu, Y. Kato, M. Ishino, M. Tanaka, N. Hasegawa, M. Nishikino, T. Ohba, T. Kaihori, Y. Ochi, T. Imazono, T. Kawachi. Soft picosecond X-ray laser nanomodification of gold and aluminum surfaces // Appl. Phys. B: Lasers and Optics, 2014, URL http://dx.doi.org/10.1007/s00340-014-5789-y. f201

[19] D.S. Ivanov, A.I. Kuznetsov, V. P. Lipp, B. Rethfeld, B.N. Chichkov, M. E. Garcia, W. Schulz. Short laser pulse nanostructuring of metals: direct comparison of molecular dynamics modeling and experiment // Appl. Phys. A, 2013. Vol. 111, p. 675-687. f201

[20] C. Wu, L.V. Zhigilei. Microscopic mechanisms of laser spallation and ablation of metal targets from large-scale molecular dynamics simulations // Appl. Phys. A, 2014. Vol. 114, p. 11-32. f201

[21] Z. Insepov, J. Rest, A.M. Yacout, A. Yu. Kuksin, G.E. Norman, V.V. Stegailov, S. V. Starikov, A. V. Yanilkin. Derivation of kinetic coefficients by atomistic methods for studying defect behavior in Mo //J. Nucl. Mat., 2012. Vol. 425, no. 1-3, p.41-47, URL http://linkinghub.elsevier.com/retrieve/ pii/S0022311511007951. f203, 205, 206

[22] M. S. Veshchunov, A. V. Boldyrev, V. D. Ozrin, V. E. Shestak, V. I. Tarasov, G.E. Norman, A. Yu. Kuksin, V.V. Pisarev, D. E. Smirnova, S.V. Starikov, V. V. Stegailov, A. V. Yanilkin. Development of the Mechanistic Fuel Performance and Safety Code SFPR Using the Multi-Scale Approach // TMS2013 Supplemental Proceedings: John Wiley & Sons, Inc., 2013, p. 655664. f203, 205, 206

[23] S.V. Starikov, Z. Insepov, J. Rest, A.Yu. Kuksin, G.E. Norman, V.V. Stegailov, A. V. Yanilkin. Radiation-induced damage and evolution of defects in Mo // Phys. Rev. B, 2011. Vol. 84, no. 10, p. 104109. f204, 205, 206

[24] D. E. Smirnova, S.V. Starikov, V.V. Stegailov. Interatomic potential for uranium in a wide range of pressures and temperatures // J. Phys. Condens. Matter, 2012. Vol. 24, no. 1, p. 015702. f204, 205

[25] D.E. Smirnova, A. Yu. Kuksin, S.V. Starikov, V.V. Stegailov, Z. Insepov, J. Rest, A.M. Yacout. A ternary EAM interatomic potential for U-Mo alloys with xenon // Model. Simul. Mater. Sci. Eng., 2013. Vol. 21, no.3, p. 035011. f204, 205, 206

[26] А. Ю. Куксин, Д. Е. Смирнова. Расчет коэффициентов диффузии дефектов и ионов в UO2 // ФТТ, 2014. Т. 56, №6, c. 1166-1175. f205, 206

[27] С. В. Стариков. Атомистическое моделирование образования дефектов при пролете осколков деления в UO2 // ТВТ, 2014 (в печати). f205, 206

[28] A.Yu. Kuksin, G.E. Norman, V.V. Pisarev, V.V. Stegailov, A.V. Yanilkin. Theory and molecular dynamics modeling of spall fracture in liquids // Phys. Rev. B, 2010. Vol. 82, no. 17, p. 174101. f207, 209, 210

[29] Z.-J. Wang, C. Valeriani, D. Frenkel. Homogeneous bubble nucleation driven by local hot spots: a molecular dynamics study //J. Phys. Chem. B, 2009. Vol. 113, no. 12, p. 3776-3784. f207

[30] J. W. P. Schmelzer. Crystal nucleation and growth in glass-forming melts: Experiment and theory // J. Non. Cryst. Solids, 2008. Vol. 354, no. 2-9, p. 269-278. f207

[31] Г. Э. Норман, В. В. Стегайлов. Гомогенная нуклеация в перегретом кристалле. Молекулярно-динамический расчет // ДАН, 2002. Т. 386, № 3, с. 328-332. f209

[32] A.Yu. Kuksin, I.V. Morozov, G.E. Norman, V.V. Stegailov, I. A. Valuev. Standards for molecular dynamics modelling and simulation of relaxation // Mol. Simul., 2005. Vol. 31, no. 14-15, p. 1005-1017. f209

[33] T. T. Bazhirov, G. E. Norman, V. V. Stegailov. Molecular dynamics simulation of cavitation in a lead melt at negative pressures // Russian Journal of Physical Chemistry, 2006. Vol. 80, no. S1, p. S90-S97. f209

[34] А. Ю. Куксин, А. В. Янилкин. Кинетическая модель разрушения при высокоскоростном растяжении на примере кристаллического алюминия // ДАН, 2007. Т. 413, № 5, с. 615-619. f210

[35] A.Yu. Kuksin, P. R. Levashov, V.V. Pisarev, M. E. Povarnitsyn, A.V. Yanilkin, A. S. Zakharenkov. Model of fracture of liquid aluminum based on atomistic simulations // Physics of Extreme States of Matter—2011, 2011, p. 57. f210

[36] D. Rossinelli, B. Hejazialhosseini, P. Hadjidoukas, C. Bekas, A. Curioni, A. Bertsch, S. Futral, S. J. Schmidt, N. A. Adams, P. Koumoutsakos. 11 PFLOP/s Simulations of Cloud Cavitation Collapse // Proceedings of SC13: International Conference for High Performance Computing, Networking, Storage and Analysis. SC '13: ACM, 2013, p. 3:1-3:13. f210

[37] В. В. Писарев. Определение свободной энергии поверхности кристалл-расплав // ТВТ, 2012. Т. 50, №6, c. 769-774. f210

[38] G. S. Smirnov, V. V. Stegailov. Melting and superheating of si methane hydrate: Molecular dynamics study // J. Chem. Phys., 2012. Vol. 136, no. 4, p. 044523. f211

[39] G.S. Smirnov, V.V. Stegailov. Toward Determination of the New Hydrogen Hydrate Clathrate Structures // J. Phys. Chem. Lett., 2013. Vol. 4, p. 35603564. f211

[40] N. English. Massively-Parallel Molecular Dynamics Simulation of Clathrate Hydrates on Blue Gene Platforms // Energies, 2013. Vol. 6, no. 6, p. 30723081. f212

[41] S.J. Stuart, A. B. Tutein, J. A. Harrison. A reactive potential for hydrocarbons with intermolecular interactions // J. Chem. Phys., 2000. Vol. 112, no. 14, p. 6472. f215

[42] S.L. Mayo, B. D. Olafson, W. A. Goddard. DREIDING: a generic force field for molecular simulations // J. Phys. Chem., 1990. Vol. 94, no. 26, p. 8897-8909. f215

[43] Н. И. Хохлов, И. Б. Петров. Решение больших задач сейсмики на высокопроизводительных вычислительных системах: Тезисы докладов: НСКФ, 2013. f216, 227

[44] В. Е. Фортов, А. Г. Храпак, С. А. Храпак, В. И. Молотков, О. Ф. Петров. Пылевая плазма // УФН, 2004. Т. 174, №5, c. 495-544. f216, 218

[45] G. Norman, V. Stegailov, A. Timofeev. Abnormal Kinetic Energy of Charged Dust Particles in Plasmas // Contrib. to Plasma Phys., 2010. Vol. 50, no. 1, p. 104-108. f217

[46] Г. Э. Норман, В. В. Стегайлов, А. В. Тимофеев. Аномальная кинетическая энергия системы пылевых частиц в плазме газового разряда // ЖЭТФ,

2011. Т. 140, №5, c. 1017-1032. f217

[47] G. E. Norman, A. V. Timofeev. Kinetic temperature of dust particle motion in gas-discharge plasma // Phys. Rev. E, 2011. Vol. 84, no. 5, p. 056401. f217

[48] G. E. Norman, A. V. Timofeev. Anomalous kinetic energy of charged dust particles in gas discharge plasma // Ukr. J. Phys., 2011. Vol. 56, no. 12, p. 1300-1304. f217

[49] Г. Э. Норман, А. В. Тимофеев. Применение понятия «температура» для описания динамики пылевых частиц в плазме газового разряда // ДАН,

2012. Т. 446, №4, c. 393-397. f217

[50] А. В. Ланкин. Столкновительная рекомбинация в неидеальной плазме // ЖЭТФ, 2008. Т. 134, c. 1013-1023. f220

[51] A. V. Lankin, G.E. Norman. Collisional recombination in strongly coupled, plasmas // J. Phys. A Math. Theor., 2009. Vol. 42, no. 21, p. 214042. f220

[52] A. Lankin, G. Norman. Density and Nonideality Effects in Plasmas // Contrib. to Plasma Phys., 2009. Vol. 49, no. 10, p. 723-731. f220

[53] A. Lankin, G. Norman. Recombination in Dense Ion Plasmas // Contrib. to Plasma Phys., 2013. Vol. 53, no. 10, p. 711-720. f220

[54] А. В. Ланкин, Г. Э. Норман, В. В. Стегайлов. Атомистическое моделирование взаимодействия электролита с графитовыми наноструктурами в перспективных суперконденсаторах // ТВТ, 2010. Т. 48, № 6, c. 877885. f223

[55] A. Bhatele, L. V. Kale, S. Kumar. Dynamic Topology Aware Load Balancing Algorithms for Molecular Dynamics Applications // Proceedings of the 23rd International Conference on Supercomputing. ICS '09: ACM, 2009, p. 110-116. f227

[56] В. В. Чуданов, С. А. Горейнов, А. Е. Аксенова, В. А. Первичко, А. А. Макаревич. Новый метод решения CFD .задач на кластерных ЭВМ петафлопсной производительности, 2013. Т. 5, № 1, c. 3-14, URL http: //psta.psiras.ru/read/psta2014_1_3-14.pdf. f227

[57] П. А. Жиляев, В. В. Стегайлов. Ab initio молекулярная динамика: перспективы использования многопроцессорных и гибридных суперЭВМ // Вычислительные методы и программирование, 2012. Т. 13, №2, c.37-45. f228

[58] I. Bethune, A. Carter, X. Guo, P. Korosoglou. Million Atom KS-DFT with CP2K // PRACE White Paper, 2011. f229

[59] O. Ayala, L.-P. Wang. Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition // Parallel Computing, 2013. Vol. 39, no. 1, p. 58-77. f230

[60] T. Hoefler, M. Snir. Generic Topology Mapping Strategies for Large-scale Parallel Architectures // Proceedings of the International Conference on Supercomputing ICS '11: ACM, 2011, p. 75-84. f231

[61] ExaScale Computing Study: Technology Challenges in Achieving Exascale Systems. DARPA Inf. Processing Techn. Office, 2008, URL http://www.cse. nd.edu/Reports/2008/TR-2008-13.pdf. f231

[62] S. Plimpton. Fast Parallel Algorithms for Short-Range Molecular Dynamics // J. Comput. Phys., 1995. Vol. 117, no. 1, p. 1-19. f232

[63] J. A. Anderson, C.D. Lorenz, A. Travesset. General purpose molecular dynamics simulations fully implemented on graphics processing units // J. Comput. Phys., 2008. Vol. 227, no. 10, p. 5342-5359. f232

[64] J.E. Stone, D.J. Hardy, I. S. Ufimtsev, K. Schulten. GPU-accelerated molecular modeling coming of age // J. Mol. Graph. Model., 2010. Vol. 29, no. 2, p. 116-125. f232

[65] M. S. Friedrichs, P. Eastman, V. Vaidyanathan, M. Houston, S. Legrand, A. L. Beberg, D. L. Ensign, C.M. Bruns, V. S. Pande. Accelerating molecular dynamic simulation on graphics processing units // J. Comput. Chem., 2009. Vol. 30, no. 6, p. 864-872. f232

[66] P. H. Colberg, F. Hofling. Highly accelerated simulations of glassy dynamics using GPUs: Caveats on limited floating-point precision // Comput. Phys. Commun., 2011. Vol. 182, no. 5, p. 1120-1129. f232

[67] I. V. Morozov, A.M. Kazennov, R. G. Bystryi, G.E. Norman, V.V. Pisarev, V. V. Stegailov. Molecular dynamics simulations of the relaxation processes in the condensed matter on GPUs // Comput. Phys. Commun., 2011. Vol. 182, no. 9, p. 1974-1978. f232

[68] R. G. Bystryi, I. V. Morozov. GPU-accelerated molecular dynamics simulations of nonideal plasmas // Physics of Extreme States of Matter-2012. --Chernogolovka: IPCP RAS, 2012, p. 134-136. f234

[69] W. M. Brown, A. Kohlmeyer, S. J. Plimpton, A. N. Tharrington. Implementing molecular dynamics on hybrid high performance computers — Particle-particle particle-mesh // Comput. Phys. Commun., 2012. Vol. 183, no. 3, p. 449-459. f234

[70] С. С. Андреев, С. А. Дбар, А. О. Лацис, Е. А. Плоткина. Макет гибридной реконфигурируемой вычислительной системы и реализация на нем вычислений с повышенной точностью: Тезисы докладов: НСКФ, 2013. f235

[71] http://top500.org/lists/2013/11/. f236

Об авторах:

Алексей Юрьевич Куксин

Кандидат физико-математических наук, старший научный сотрудник ОИВТ РАН, доцент кафедры молекулярной физики МФТИ, лауреат гранта Президента РФ для молодых учёных-кандидатов наук 2011 г. e-mail: alexey.kuksin@gmail.com

Александр Валерьевич Ланкин

Кандидат физико-математических наук, старший научный сотрудник ОИВТ РАН, ассистент кафедры общей физики МФТИ, лауреат медали РАН с премией для студентов вузов России 2007 г.

e-mail: alex198508@yandex.ru

Игорь Владимирович Морозов

Кандидат физико-математических наук, заведующий лабораторией ОИВТ РАН, доцент кафедр информатики МФТИ и прикладной математики ВШЭ, лауреат медали РАН с премией для молодых учёных России 2004 г. и гранта Президента РФ для молодых учёных-кандидатов наук 2010 г.

e-mail: morozov@ihed.ras.ru

Генри Эдгарович Норман

Доктор физико-математических наук, главный научный сотрудник ОИВТ РАН, профессор кафедры молекулярной физики МФТИ, Соросовский профессор 1997 и 1998 г., гранты Москвы «Профессор—2002, 2004, 2005» e-mail: genri.norman@gmail.com

Никита Дмитриевич Орехов

Стажёр-исследователь ОИВТ РАН, студент МФТИ e-mail: nikita.orekhov@gmail.com

Василий Вячеславович Писарев

Кандидат физико-математических наук, старший научный сотрудник ОИВТ РАН, ассистент кафедры молекулярной физики МФТИ, лауреат медали РАН с премией для студентов вузов России 2011 г.

e-mail: pisarevvv@gmail.com

Григорий Сергеевич Смирнов

Младший научный сотрудник ОИВТ РАН, аспирант МФТИ, лауреат медали РАН с премией для студентов вузов России 2013 г.

e-mail: grsmirnov@gmail.com

Сергей Валерьевич Стариков

Кандидат физико-математических наук, заведующий лабораторией ОИВТ РАН, ассистент кафедры общей физики МФТИ, лауреат гранта Президента РФ для молодых учёных-кандидатов наук 2012 г.

e-mail: redshuhart@yandex.ru

Владимир Владимирович Стегайлов

Доктор физико-математических наук, заведующий отделом ОИВТ РАН, доцент кафедры молекулярной физики МФТИ, лауреат медали РАН с премией для студентов вузов России 2004 г., лауреат гранта Президента РФ для молодых учёных-кандидатов наук 2007 и 2010 гг. e-mail: stegailov@gmail.com

Алексей Владимирович Тимофеев

Кандидат физико-математических наук, заведующий лабораторией ОИВТ РАН, доцент кафедры прикладной математики ВШЭ, ассистент кафедры общей физики МФТИ, лауреат медали РАН с премией для студентов вузов России 2010 г.

^ e-mail: timofeevalvl@gmail.com

Образец ссылки на эту публикацию:

А. Ю. Куксин, А. В. Ланкин, И. В. Морозов, Г. Э. Норман, Н. Д. Орехов, В. В. Писарев, Г. С. Смирнов, С. В. Стариков, В. В. Стегайлов, А. В. Тимофеев. ЗАЧЕМ и КАКИЕ нужны суперкомпьютеры экс-афлопсного класса? Предсказательное моделирование свойств и многомасштабных процессов в материаловедении // Программные системы: теория и приложения: электрон. научн. журн. 2014. T. 5, № 1(19), с. 191-244.

URL: http://psta.psiras.ru/read/psta2014_1_191-244.pdf

A. Y. Kuksin, A. V. Lankin, I. V. Morozov, G. E. Norman, N. D. Orekhov, V. V. Pisarev, G. S. Smirnov, S. V. Starikov, V. V. Stegailov, A. V. Timofeev. Predictive modeling and simulation of properties and multi-scale processes in materials science. Tasks for Exaflops-era supercomputers.

Abstract. The approach is developed which allows to find out the problems which need for their solution exaflops supercomputers. The approach is demonstrated at the examples of topical problems of material science, condense matter and dense plasma physics where atomistic modeling is necessary to apply. The correspondence is established for each problem between phenomena studied and computational cores number needed. Modeling parallel programs scalability is shown as well as perspective of the modeling methods predictive ability extension with the increase of computational cores number and / or use of special architecture (graphical processing units).

The following problems are considered: 1) surface modification at processing of metals by sub-picosecond laser pulses, 2) radiation-induced aging of nuclear reactors fuels, 3) phase transition kinetics in metastable liquids, 4) methane and hydrogen gas hydrates structures and computation of their properties, 5) polymers multiscale models, 6) dusty plasmas, 7)ion recombination in liquid and gaseous dielectric media at discharge break and relaxation, 8) electric double layer between graphite and electrolyte, influence of electron-hole electrode structure on capacity. Predictive modeling reliability is checked by comparisons with experiments.

The modeling methods hierarchy, which is necessary to describe properties of matter at different space and time scales, is considered in frames of the multiscale approach. Density functional theory (quantum molecular dynamics) is applied at the deepest nm/pm scale to model electron dynamics and to construct effective interaction potentials between particles. Classical molecular dynamics modeling is used to treat moving atoms systems up to micro-scale. Kinetic theory and continuum mechanics is used to proceed with micro-scale. Particular attention is paid to the exchange of information between different scales, i.e. to the unified description of systems from nano to micro levels. Parallelization efficiency comparison is performed for three classes of problems at fat tree and torus topologies (in Russian).

Key Words and Phrases: atomistic modeling, electronic structure, molecular dynamics, multiscale modelling, radiation aging, laser ablation, nucleation, hydrates, polymers, dusty plasma, electrochemistry, parallel efficiency.