Содержание к диссертации
Введение
Глава 1. Обзор литературы 8
1.1. Метод молекулярной динамики 8
1.2. Структурные свойства газовых гидратов 17
1.3. Направления исследований газовых гидратов 19
Глава 2. Гидраты метана 24
2.1. Молекулярно-динамическая модель 24
2.2. Расчет кривой плавления структуры I гидрата метана 28
2.3. Распад перегретой структуры I 32
2.4. Выводы 38
Глава 3. Гидраты водорода 39
3.1. Введение 39
3.2. Устойчивость структур 44
3.3. Диффузия молекул водорода 52
3.4. Анализ средних квадратов смещений молекул водорода 55
3.5. Анализ траекторий отдельных молекул 57
3.6. Выводы 63
Глава 4. Применение суперкомпьютеров для моделирования водных систем 64
4.1. Введение 64
4.2. Сравнение современных ускорителей 67
4.3. Реализация молекулярно-динамических алгоритмов на ускорителях 73
4.4. Результаты тестов
4.5. Выводы 81
Основные результаты и выводы 83
Список сокращений и условных обозначений 84
Литература
- Структурные свойства газовых гидратов
- Расчет кривой плавления структуры I гидрата метана
- Диффузия молекул водорода
- Сравнение современных ускорителей
Введение к работе
Актуальность темы исследования. Гидраты природных газов, или газовые гидраты — это нестехиометрические соединения включения, в которых молекулы газа (молекулы-гости) заключены в полостях трёхмерной решётки из молекул воды (каркас хозяина). Стабилизация водных клатратных каркасов, термодинамически менее стабильных, чем лёд или жидкая вода при тех же условиях, обеспечивается за счёт взаимодействий гость-хозяин.
Для извлечения углеводородов из газогидратных пластов необходимо знание их физических характеристик, таких как фазовая диаграмма, кинетика разложения, фильтрационные свойства и др. Существующие математические модели хорошо описывают фазовые диаграммы и свойства гидратов метана и других простых углеводородов в широком диапазоне давлений и температур. Однако на зону стабильности гидратов сильно влияет наличие примесей других газов, пористость и проницаемость газосодержащих пород, переход вода-лед при отрицательных температурах. На кривую фазового равновесия также влияет наличие примесных солей в воде, что критично при разработке морских месторождений гидратов. Указанные эффекты могут быть учтены при разработке многомасштабных моделей формирования и разрушения гидратов.
Хранение и транспортировка водорода в гидратах также требует детального знания фазовой диаграммы данной системы. Молекулы водорода из-за малого размера при высоких давлениях легко растворяются в различных формах водного льда и клатратных структур. Переходы между такими структурами могут сопровождаться изменением равновесной концентрации молекул водорода, что нежелательно с практической точки зрения.
Цели и задачи исследования:
-
Сравнение точности различных потенциалов взаимодействия для моделирования гидратов метана и водорода.
-
Расчет кривой плавления структуры I гидрата метана при давлениях
выше 40 МПа. Определение ее кинетической границы устойчивости, в том числе при варьировании скорости нагрева и степени заполнения полостей.
-
Анализ устойчивости предполагаемых на основе экспериментальных данных вариантов структуры новой фазы гидрата водорода в предполагаемой области ее существования.
-
Определение характера диффузии газа в Со и sT' структурах гидрата водорода.
-
Оценка максимально возможного времени расчетов классических моле-кулярно-динамических моделей водных систем с использованием суперкомпьютеров и сравнение эффективности соответствующих высокопроизводительных вычислительных систем.
Научная новизна. В работе впервые рассчитана кривая плавления гидрата метана в диапазоне давлений от 0 до 0,5 ГПа с использованием моделей SPC/E, TIP4P/2005, Т1Р4Р/1се. Получена зависимость кинетической границы устойчивости перегретого гидрата метана от скорости нагрева и степени заполнения полостей. Сделана проверка устойчивости на наносекундных временах структур а-кварца, Со, sT' в диапазоне давлений 0,5-0,7 ГПа и температур 170-250 К для сочетания моделей TIP4P/2005, Т1Р4Р/1се для воды и двух моделей для водорода. Произведена оценка систематической погрешности определения области существования устойчивых структур на фазовой диаграмме. Представлены результаты расчета диффузии молекул-гостей в гидратах водорода на наносекундных временах, количественно характеризующие ее анизотропный и аномальный характер. Впервые введена метрика для ранжирования высокопроизводительных вычислительных систем «полная пиковая вычислительная производительность - время расчета» и ее использование для сравнения эффективности решения на современных суперкомпьютерах типичных классических молекулярно-динамических моделей водных систем.
Теоретическая и практическая значимость работы. Газовые гидраты рассматриваются как потенциальный источник углеводородов, а также
как средство хранения и транспортировки газов. Данные по кривой плавления гидратов метана позволяют оценить точность классических потенциалов для воды при моделировании различных свойств газовых гидратов. Исследования структур гидратов водорода позволили уточнить экспериментальные данные о фазовой диаграмме данной системы. Сравнение эффективности суперкомпьютеров определяет диапазон явлений, доступных для моделирования.
Методы исследования. Молекулярная динамика применяется для моделирования газовых гидратов более тридцати лет. В основе метода лежит численное интегрирование классических уравнений движения Ньютона для системы N частиц. Силы между атомами определяются из эмпирических потенциалов взаимодействия, которые позволяют качественно и количественно воспроизводить свойства газовых гидратов.
Положения, выносимые на защиту:
-
Результаты расчета температур плавления гидрата метана в диапазоне давлений от 0 до 0,5 ГПа с использованием моделей SPC/E, ТІР4Р/2005, ТІР4Р/ Ice. Зависимость кинетической границы устойчивости перегретого гидрата метана от скорости нагрева и степени заполнения полостей.
-
Проверка устойчивости на наносекундных временах структура-кварца, Со, sT' гидратов водорода в диапазоне давлений 0,5-0,7 ГПа и температур 170-250 К для сочетания моделей TIP4P/2005, Т1Р4Р/1се для воды и двух моделей для водорода. Оценка систематической погрешности определения области существования устойчивых структур на фазовой диаграмме.
-
Результаты расчета среднеквадратичных смещений молекул-гостей в гидратах водорода на наносекундных временах, показывающие анизотропный и аномальный характер диффузии.
-
Метрика для ранжирования высокопроизводительных вычислительных систем «полная пиковая вычислительная производительность - время расчета» и ее использование для сравнения эффективности решения на современных суперкомпьютерах классических молекулярно-динамических моделей водных
систем.
Степень достоверности полученных результатов. Результаты моле-кулярно-динамичеких расчетов находятся в согласии с известными экспериментальными данными, а также расчетными данными других авторов.
Апробация результатов. Основные результаты диссертации докладывались на конференциях: «Nucleation Theory and Applications» (Дубна, Россия, 2010, 2013, 2015), «Interaction of Intense Energy Fluxes with Matter» (Эльбрус, Россия, 2011, 2013), научные конференции МФТИ (Москва, Россия, 2009-2014), «Проблемы физики ультракоротких процессов в сильнонеравновесных средах» (Новый Афон, Абхазия, 2010, 2011, 2012), 7th International Conference on Gas Hydrates (Эдинбург, Великобритания, 2011), 8th International Workshop on Methane Hydrate Research Sz Development (Саппоро, Япония, 2012), CPMD-Meeting (Лейпциг, Германия, 2013), EMLG-JMLG annual meeting (Лилль, Франция, 2013), Water Europe (Сарагоса, Испания, 2014), Dushanbe Symposium on Computational Materials and Biological Sciences (Душанбе, Таджикистан, 2014), XIV российская конференция по теплофизическим свойствам веществ (Казань, Россия, 2014), International Supercomputing Conference (Франкфурт, Германия, 2015), International Conference on Computer Simulation in Physics and beyond (Москва, Россия, 2015), «Суперкомпьютерные дни в России» (Москва, Россия, 2015).
Публикации. По материалам диссертации опубликовано 9 печатных работ. Работ, опубликованных в рецензируемых научных изданиях, рекомендованных ВАК — 5.
Личный вклад автора. Все представленные в диссертации результаты получены лично автором. Задачи численных экспериментов по диссертационной работе сформулированы под руководством В. В. Стегайлова. Выводы и основные положения, выносимые на защиту, сформулированы автором.
Структура и объем диссертации. Диссертация состоит из введения, четырех глав, заключения, списка сокращений и библиографии. Общий объем диссертации 104 страницы, включая 27 рисунков. Библиография включает 165
5 наименований.
Структурные свойства газовых гидратов
Большая часть вычислительного времени в молекулярной динамике тратится на расчет сил. Их эффективное вычисление для больших систем требует использования специальных алгоритмов [9].
Для расчета короткодействующих сил используется два основных метода. Первый метод («списки Верле») заключается в составлении списков соседей ближайших атомов. Обычно при составлении списков сохраняются все ближайшие атомы на расстоянии ra = rc + S. Список используется несколько шагов интегрирования, далее он составляется снова, чтобы атомы на расстоянии г rs не могли попасть в область г гс. Выбор параметра 5 зависит от конкретных условий расчета (температура, плотность и т.п.), хотя он должен быть много меньше расстояния гс.
Второй метод заключается в разбиении расчетной ячейки на трехмерные области размером d гс. Это позволяет искать ближайших соседей атома только в 27 областях, что уменьшает время на построение списков Верле.
При использовании нескольких процессоров требуется параллелизация алгоритма. Наиболее простой для реализации способ — декомпозиция по атомам. В этом случае каждый из Р процессоров в начале расчета получает группу N/ Р атомов. Процессор вычисляет силы, а затем новые координаты и скорости только этих атомов, независимо от их положения в пространстве. Метод прост в реализации, однако требует глобальных обменов данными между процессорами, поэтому используется преимущественно на системах с общей памятью. Дополнительный выигрыш в скорости может быть получен при использовании третьего закона Ньютона F = — Fjj. В этом случае при масштабировании числа атомов и/или числа процессоров время на вычисления растет как тур + N, время на обмен данными и объем памяти — как 2N.
Более сложный подход заключается в пространственной декомпозиции атомов. Расчетная ячейка разбивается на подъячейки, по одной на процессор. Про 15 цессор вычисляет силы, координаты и скорости всех атомов в данной подъячей-ке. При переходе из одной области в другую атом переназначается другому процессору. Размер и форма областей зависит от количества атомов, числа процессоров и размеров расчетной ячейки. Они должны подбираться таким образом, чтобы каждая область имела форму, близкую к кубической. Это позволяет снизить обмен данными между процессорами. Алгоритм масштабируется нели N , г /W\2/3 г неино: время на вычисления растет как 2гр + brs (-р-) , на обмен данными как 6rs ( ) , затраты на память как +6rs ( ) . В то же время, при неравномерной плотности атомов в расчетной ячейке бывает выгоднее с вычислительной точки зрения выполнять неравномерное разбиение ячейки, чтобы каждый процессор выполнял вычисления для примерно одинакового количества атомов.
Пространственное разбиение используется и для распараллеливания алгоритма РРРМ [10]. Каждый процессор обрабатывает точки сетки, лежащие в его пространственной области. Дополнительно к ним включаются фиктивные ближайшие точки соседних областей для корректного учета зарядов на границе. На первом шаге каждый процессор выполняет интерполяцию зарядов на сетку и производится суммирование зарядов с фиктивных точек.
Для наиболее эффективного вычисления трехмерного преобразования Фурье каждый процессор должен последовательно обрабатывать подмножество одномерных столбцов трехмерной сетки. Например, проекция сетки на плоскость ху делится на квадраты по числу процессоров, каждый процессор обрабатывает все точки по направлению z для своего квадрата. Это требует дополнительного обмена данными, но позволяет ускорить расчет. Таким образом, на втором шаге решается уравнения Пуассона с помощью быстрого преобразования Фурье и учетом вышеуказанного замечания. Далее производится расчет поля на сетке и расчет сил на заряды при помощи интерполяции. 1.1.5. Расчёт термодинамических величин
При аналитическом решении дифференциальных уравнений движения при заданных начальных условиях справедлива теорема существования и единственности решения задачи Коши. Однако численные методы ведут к возникновению неустойчивого решения [11], то есть решение задачи Коши с близкими начальными условиями дает экспоненциально разбегающиеся траектории: \Ar{t)\ exp{\t). (1.20)
Ньютоновская динамика сохраняется только на временах, меньших времени динамической памяти tdm- Оно определяет промежуток времени, за который теряется корреляция численного и точного решений для одинаковой начальной конфигурации. Таким образом, точное предсказание траектории атомов возможно только в пикосекундном интервале и имеет смысл говорить только о средних (по ансамблю или времени) величинах. Основным является микроканонический ансамбль NVE, в котором сохраняется число частиц N, объем системы V и полная энергия (Е), равная сумме кинетической (К) и потенциальной энергии (U) атомов. Из-за использования конечно-разностной аппроксимации значение полной энергии не сохраняется, а флуктуирует около своего среднего значения, поэтому ансамбль NVE в молекулярной динамике не полностью эквивалентен термодинамическому ансамблю NVE.
Расчет кривой плавления структуры I гидрата метана
В данной работе исследовалась граница трёхфазного равновесия в гидратах метана (кристалл-жидкость-газ) [65, 66]. Аналогичные вычисления были сделаны в работе Конде и Беги [67], где проводилось длительные (до 1 мкс) вычисления молекулярно-динамических траекторий в NPT ансамбле. В начале расчёта в одной половине прямоугольной расчётной ячейке находится кристалл гидрата метана, в другой — его расплав. Рассмотрим процессы, которые могут происходить в такой системе в ходе молекулярно-динамического расчёта. Поскольку точка фазового равновесия заранее неизвестна, в системе наблюдается рост кристалла, либо его плавление. Такой процесс зависит от кинетической энергии всех молекул. Поскольку полная энергия не является постоянной величиной, потенциальная энергия изменяется при фазовом переходе. Уменьшение потенциальной энергии соответствует росту кристалла, а её рост — плавлению. Проводя серию расчётов с различной температурой, но постоянным давлением, можно определить равновесную температуру как среднее между наибольшей Т, при которой наблюдается рост кристалла, и наименьшей Т, при которой наблюдается плавление. Главным недостатком этого метода является необходимость проведения большого количества длительных расчётов. По этой причине в данной работе использовалась другая методика [68, 69].
В качестве начального состояния рассматривалась система, состоящая из п х п х 2п элементарных ячеек, где п менялось от 2 до 6 (подробное обсуждение размерного эффекта приведено ниже). Одна половина расчётной ячейки плавилась, в то время как другая половина оставалась «замороженной». После этого в ходе короткого NPT-расчёта система выводилась в желаемый диапазон температур и давлений. Далее проводился молекулярно-динамический расчёт в NVE-ансамбле, в ходе которого система выходила на равновесие. В зависимости от начальной кинетической энергии возможны два процесса: рост кристалла, либо его плавление. Если начальная температура системы слишком высокая (или, наоборот, низкая) кристалл полностью расплавится (или, наоборот, вся жидкость кристаллизуется). Такой процесс не позволит определить точку фазового равновесия, поскольку в процессе расчёта должна образоваться равновесная трёхфазная система кристалл -жидкость- газ. Именно поэтому степени перегрева или переохлаждения не должны быть слишком большими. Пример трёхфазной системы, полученной в ходе моделирования, показан на рисунке 2.2. Когда объем фазы I перестает изменяться, температура и давление соответствуют равновесным. Пример изменения давления и температуры в ходе расчета приведен на рисунке 2.3.
Анализ показывает, что молекулы метана собираются в пузырь цилиндри зо Время, не Рис. 2.3. Пример эволюции температуры (нижняя кривая) и давления (верхняя кривая) в ходе установления фазового равновесия. Горизонтальные линии показывают равновесные значения для конкретной МД траектории. ческой формы. Именно по этой причине изучение роста газовых гидратов сильно затруднено — слишком мало молекул метана находится вблизи границы раздела кристалл-жидкость. По оценке Конде и Беги, скорость роста кристалла о составляет 0,8 А/нс. Авторами исследовались только перегретые состояния, что позволяет добиться уменьшения времени расчёта.
Типичная расчётная ячейка в предыдущих работах других авторов содержит 2x2x4 элементарных ячеек. Известно, что размерные эффекты малы для простых кристаллов, таких как алюминий или лёд. Для газовых гидратов такой анализ раньше не проводился. Поэтому была исследована зависимость температуры плавления от размера системы (п х п х 2 п, где п менялось от 3 до 6) для Т1Р4Р/1се модели.
Численные значения для Т1Р4Р/1се модели и различных размеров систем приведены в таблице 2.4. Сходимость достигается при использовании п = 5 элементарных ячеек. Аналогичные результаты были получены и для других моделей воды.
Фазовая диаграмма вода-метан. Чернвіе треуголвники — эксперименталвная кривая трехфазного равновесия кристалл-жидкоств-газ гидрата метана [70]. Синие звездві — резулвтатві расчёта Йенсена и др. [71] для ТІР4Р/ІСЄ модели. Заполненнвіе голубвіе, зеле-нвіе, краснвіе и оранжеввіе символві показвівают резулвтатві из работ [67] и [72]. Пуствіе символві — резулвтатві данной работві. Оранжеввіе треуголвники — моделв SPC/E, зеленвіе ромбві и краснвіе квадратві — ТІР4Р/2005 с % = 1, 07 и 1,00 соответственно, голубвіе круги — ТІР4Р/Ісе.
По данным Конде и Веги [67], модель ТІР4Р/Ісе дает наилучшее согласие с экспериментальными данными. Аналогичные вычисления делались в работе Йенсена и др.[71] путем вычисления свободной энергии методом Монте-Карло. Их результаты дают худшее согласие с экспериментом (хотя потенциал для метана несколько отличался). Результат для Т1Р4Р/1се модели, полученный в данной работе, согласуется с данными Йенсена, хотя данные для ТІР4Р/2005 соответствуют работе Конде и Веги. Возможно, это связано с меньшими размерами поперечного сечения в работе Конде и Веги. Результаты данной работы показы
Диффузия молекул водорода
Методами молекулярной динамики были исследованы [125] предложенные экспериментаторами структуры. Структура а-кварца содержит 9 молекул воды в элементарной ячейке (позиции За, х = 0,4699 и 6с с х = 0,4141, у = 0,2681, z = 0,1188). При использовании параметров решетки, предложенных в рабо о о те [100], а = 6,24 А и с = 6,18 А среднее расстояние между атомами кислорода меньше 2 А. Это приводит к возникновению нереалистично большого давления и практически мгновеному распаду структуры в ходе расчета. Увеличение параметров а и с приведет к расхождению с экспериментальной дифракционной картиной. Также оказался нестабилен вариант структуры Q) с молекулами воды в каналах, что подтверждает предположение [101]. В то же время структуры Со (только с водородом в каналах) и sT (приведены на рис. 3.2) остаются стабильными в ходе молекулярно-динамических расчетов.
В общем случае наиболее устойчивая при данных условиях фаза имеет наименьшую свободную энергию. Однако вычисление свободной энергии для многокомпонентной системы представляет собой довольно трудную задачу [71]. В потенциальную энергию входит взаимодействие вода-вода, вода-водород и водород-водород. Структуры содержат различное число молекул на элементарную ячейку, поэтому необходимо сравнивать их энергии на 1 атом. Для первоначальной оценки была рассчитана потенциальная энергия соответствующих пустых решеток при нулевой температуре. Зависимость энергии на одну молекулу воды от плотности фазы приведена на рис. 3.3. Вариант фазы Q) с молекулами воды в каналах обозначен как CQ-I, без — Со-П. Потенциальные энергии структур Со-П и sT довольно близки. Структуры CQ-I И а-кварца имеют более высокую энергию, что объясняет их неустойчивость в МД-расчете.
Далее устойчивость данных структур проверялась при различных температурах и давлениях. В структуре Со-П было 1080 молекул Н2О и 540 молекул Н2, в структуре sT 900 Н2О и 300 Н2. Для этого последовательно запускались Рис. 3.2. Шаростержневые модели клатратных структур Со и sT . Красные стержни изображают водородные связи между молекулами воды. Молекулы Н2 изображены светло-коричневыми шарами. короткий 100 пс NPT и длинный 50 не NVE расчеты. Устойчивость проверялась вычислением структурных параметров порядка F3, ц, F4r [99]. Определим эти параметры для систем, содержащих воду. Пусть А — атом кислорода в молекуле воды, В и Bj — другие атомы кислорода в первой координационной сфере. Параметр порядка F3 показывает отклонение от тетраэдрической конфигурации водородных связей. Пусть угол в — угол в тройке B -A-Bj атомов кислорода в молекулах воды, тогда параметр F3 определяется как F3 = (cos(6)\cos(6)\ + со52(109,47)), (3.1) где угловые скобки означают усреднение по всем уникальным парам Д и Bj и О с; о -20 а-кварц С0 (с водой в каналах) С0 sT LU -40 -60 -80 0.4 0.8 1.2 Плотность, г/см3 1.6 Рис. 3.3. Зависимость потенциальной энергии на одну молекулу пустой решетки от плотности фазы при нулевой температуре (модель воды Т1Р4Р/1се). центральным атомам А. Параметр F3 близок к 0 для жидкой воды с тетраэдри-ческой конфигурацией и отличен от нуля в противном случае. Параметры ц и F T разработаны для анализа торсионных углов в системе водородных связей и рассчитываются по следующим формулам: ц = (cos (30)), (3.2) FAT = ((HmOi OiOj OjHn)). (3.3)
Здесь угол ф соответствует двугранному углу, образованному четверками атомов Hm-Oj. . .Oj-Hn. Использование этих параметров позволяет эффективно различать жидкое и твердое состояние, а также различать различные структуры газовых гидратов. На рисунках 3.4 и 3.5 приведены примеры зависимости указанных параметров порядка от времени для различных температур и степеней заполнения. Давление в указанных расчетах составляет 0,6 ГПа. Для Q) структуры показана зависимость для нулевой (пустая решетка) и 100% степени заполнения полостей (слева и справа соответственно). При нулевой степени заполнения и температуре 220 К решетка распадается за времена порядка 0,5 наносекунд. Это отражено в первых трех графиках рисунка 3.4. При 100% заполнении полостей решетка структуры остается стабильной в исследуемом диапазоне температур.
На рисунке 3.5 приведены зависимости указанных параметров порядка для 100 и 75% заполнения полостей в гидрате sТ при различных температурах. Видно, что хотя значения параметров F3 практически совпадают со структурой Со, различия в значениях параметров F , F41- позволяют разделять данные структуры. В рассматриваемом диапазоне температур в структуре sT фазовые превращения начинаются уже при 75% заполнении полостей и температуре 260 К — начальная структура оказывается метастабильной.
Эволюция параметров порядка F3, F4 A, F4r в фазе sT Процессы распада в гидратах СО2 и СН4 изучались в различных работах [66, 126-129]. Среднее время жизни зависит от многих факторов: степени метастабильности, числа молекул, потенциала взаимодействия. Однако даже ограниченная устойчивость структур (например, в пределах 50 не) может служить доказательством, что они могут существовать в данной области фазовой диаграммы.
Гидрат sT не распадается в ходе 50 не расчета при (140 К; 0,2 ГПа), (140; 1), (220; 0,6) и (300; 0,2). Это применительно для всех потенциалов взаимодействия, использовавшихся в расчетах. В точке (300 К; 1 ГПа) поведение системы отличается в зависимости от используемого потенциала взаимодействия. Она стабильна только при использовании TIP4F/Ice модели воды и трехточечной модели молекулы водорода, полностью распадается в конце NPT расчета для TIP4P/2005 модели и частично в ходе NVE расчета для Т1Р4Р/1се модели и UA LJ модели водорода. В последнем случае давление и температура изменяются в ходе NVE расчета из-за фазового перехода, в ходе которого часть кристалла плавится. В этом случае напряжения становятся негидростатическими, однако грубо позволяет оценить температуру плавления sT структуры. Она достаточно хорошо согласуется с кривой плавления гидратов водорода (рис. 3.6).
Гидрат Со нестабилен только в точке (300 К; 0,2 ГПа) при использовании модели воды TIP4P/2005. В 4 других точках (Т,Р) и других комбинациях силовых полей он устойчив в ходе всех 50 нс МД-расчетов.
Возможность формирования Со фазы также подтверждается ab initio расчетами в более поздней работе других авторов [134] (при этом, однако, структура sT не обсуждается). Данные результаты показывают возможность существования обоих структур гидрата водорода, хотя бы метастабильных.
Точность результатов атомистического моделирования в первую очередь зависит от точности используемых потенциалов межатомного взаимодействия. Известно, что уравнения состояния и фазовая диаграмма благородных газов хорошо описывается потенциалом Леннард-Джонса. Фазовую диаграмму воды можно качественно описать правильно выбранными потенциалами, например, TIP4P/2005 или Т1Р4Р/1се.
Сравнение современных ускорителей
Гибридные вычислительные системы, использующие несколько различных вычислительных устройств, получили широкое распространение в последние годы. Анализ списка Топ-500 лучших суперкомпьютеров мира, показывает, что наибольшей популярностью на сегодняшний день пользуются машины, использующие в качестве укорителей видеокарты NVIDIA архитектуры Kepler или сопроцессоры Intel Xeon Phi архитектуры MIC (Many Integrated Core). Несколько суперкомпьютеров используют одновременно обе архитектуры. В качестве особого примера гибридной системы можно указать суперкомпьютер с японскими 1024-ядерными процессорами PEZY-SC. Всего в ноябрьской редакции 2014 года списка Топ-500 75 машин с гибридной архитектурой. На рисунке 4.1 приведены данные по числу гибридных вычислительных систем и их пиковой производительности в разные годы.
Резкий рост количества машин с гибридной архитектурой связан в первую очередь с перспективой существенного увеличения производительности суперкомпьютера без значительного увеличения затрат. Пиковая производительность современных ускорителей достигает нескольких ТФлопс (10 операций с плавающей точкой в секунду). Однако использование значительной доли пиковой вычислительной мощности Rpeak в реальных приложениях не всегда возможно. Это связано с необходимостью полного переписывания исходного кода для обычных CPU, решения проблем с доступом к памяти ускорителя, шириной канала передачи данных между центральным процессором и ускорителем, а также особенностями работы арифметико-логических устройств ускорителей.
Предпринимаются попытки создания специализированных суперкомпьютеров для конкретных приложений, однако они не входят в список Топ-500. Для задач молекулярной динамики можно выделить суперкомпьютеры на программируемых пользователем вентильных матрицах (разновидность ПЛИСов), MDGRAPE [147] и ANTON [148] на интегральных схемах специального назначения. Использование специальной архитектуры позволяет добиться наиболее эффективного использования вычислительной мощности для конкретной задачи.
При этом «классические» многопроцессорные суперкомпьютеры, использующие для вычислений только центральные процессоры, продолжают доминировать в списке Топ-500, и на такие системы продолжают ориентироваться при разработке ПО для классической молекулярной динамики. На них проводятся рекордные по числу частиц молекулярно-динамические расчеты [149]. В сложившейся ситуации при наличии большого разнообразия аппаратного обеспечения для высокопроизводительных суперкомпьютерных расчетов особенную актуальность приобретает вопрос сравнения между собой различных альтернативных решений. Целью анализа должна быть эффективность связи 1) аппаратного обеспечения, 2) конкретной математической модели или класса моделей и 3) численных алгоритмов, в том числе с учетом уже существующего программного обеспечения (ПО) и сложности их адаптации на новые типы аппаратного обеспечения [150, 151].
Популярный тест LINPACK, использующийся для ранжирования суперкомпьютеров, не отражает особенностей многих современных алгоритмов, в том числе и алгоритмов молекулярной динамики. Для данного класса задач это приводит к существенному отличию между реальной и теоретически максимальной производительностью [152-154]. Отсюда следует, что необходим анализ эффективности работы конкретных программ и алгоритмов на различных архитектурах. В качестве примеров аппаратного обеспечения рассматриваются лучшие суперкомпьютерные системы России. Анализируемые математические модели соответствуют двум стандартным примерам: примеру леннард-джонсов-ской жидкости (типичной системе для задач статистической физики, физики конденсированных сред и физической химии) и примеру белковой молекулы в водном растворе. В качестве примера ПО мы используем пакет LAMMPS (одну из наиболее популярных и многофункциональных программных реализаций методов молекулярно-динамического моделирования) и пакеты, ориентированные на биомолекулярные исследования (NAMD и др.)
На сегодняшний день атомистическое моделирование — одно из важнейших приложений суперкомпьютерных вычислений [153]. На его основе получили развитие многомасштабные модели в химии, физике и других областях. Однако даже рекордные расчеты триллионов атомов [149] соответствуют объему всего в несколько кубических микрометров. Увеличение временных масштабов в молекулярно-динамических расчетах — еще более трудоемкая задача. Развитие таких вычислительных методов происходит одновременно с развитием суперкомпьютеров.