Доклады Российской академии наук. Химия, науки о материалах , 2022, T. 504, № 1, стр. 72-102

III.1. Межатомные потенциалы машинного обучения

Межатомные потенциалы машинного обучения (machine learning interatomic potentials, MLIP) за последнее десятилетие достигли определенной зрелости и стали достойной альтернативой обычным межатомным потенциалам. Этому успеху во многом способствовала американская Инициатива по геному материалов [38], а также поддержка Европейским советом по моделированию материалов [39], в рамках которых основаны, например, открытая и бесплатная интернет-платформа NanoHUB, база квантовых данных о свойствах материалов (OQMD) [40], автоматическая система поиска материалов (AFLOW) [41, 42] и многие другие.

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

Рис. 3.

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

Отметим, что к настоящему моменту развито уже впечатляющее количество методов машинного обучения потенциалов для атомистического моделирования и этот процесс продолжается. Причины такого состоят в том, что существует определенное противоречие в двух основных пунктах схемы на рис. 3. Если во главу угла ставить быстродействие и качество обучающей модели, то это обычно требует максимального расширения пространства дескрипторов. Такой подход можно условно назвать “математическим”. Противоположный подход, “физический”, стремится выбрать пространство дескрипторов максимально приближенным к физике процессов межчастичного взаимодействия. При этом пространство дескрипторов сужается, что накладывает ограничения на методику обучения. Найти оптимальное соотношение “математического” и “физического” подходов для различных типов материалов довольно сложно. По этой причине нам представляется важным ознакомить читателя не с описанием более чем десятка существующих методик построения MLIP, а объяснить основные принципы их работы.

Как уже отмечалось, аппроксимация вида (прогнозирование) ППЭ на основе данных первопринципных расчетов, относится к методам машинного обучения с учителем ( так называемое контролируемое обучение). Все такие методы удобно рассматривать на основе формального статистического подхода [44, 45]. Пусть у нас есть обучающий набор D из n результатов первопринципного моделирования, D = {x_i, y_i} i = 1, …, n}, где х – входной вектор параметров исследуемого материала размерности D, а у – скалярный результат расчета потенциальной энергии. Обозначим через Х матрицу размером D × n, в которой собраны все входные данные обучающей выборки, а через y – вектор выходных данных тестовой выборки. Общий результат прогнозирования можно математически сформулировать как построение некоторой функции  f, которая отображает набор обучающих данных x (серые точки на рис. 4) на прогнозируемый результат y. Наша цель – использовать такие обучающие данные, а также любые другие предшествующие знания, чтобы идентифицировать функцию, которая способна точно предсказать выходное значение для любых значений х из набора входных данных, не совпадающих с обучающим. Если допустимые выходные значения y образуют непрерывный диапазон, тогда процесс поиска функции f называется регрессией. Назовем пространством гипотез все возможные функции, которые нам может выдать процесс обучения. Пространство гипотез содержит все гипотезы (т.е. функции), которые могут быть возвращены алгоритмом обучения. Однако в этом пространстве будет много функций, которые точно воспроизводят данные обучения (рис. 4), и необходимо каким-то образом ограничить (отфильтровать) это пространство гипотез, приписывая функциям некоторые веса так, чтобы одни функции были более предпочтительными, чем другие. Например, очень часто пространство гипотез ограничивают так, чтобы оставались лишь функции, выражающие линейную зависимость (линейная регрессия) между входными и выходными значениями (сплошная линия на рис. 4). Заметим, что в пространстве гипотез может не оказаться функции, которая производит наблюдаемые выходные значения для всех возможных наборов входных значений. Это может произойти по нескольким причинам:

Рис. 4.

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

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

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

3) некоторые входные данные, имеющие особую важность для вычисления правильных выходных значений, отсутствуют (например, точка вблизи х = 0, значение которой (y(0)) могло бы полностью отвергнуть гипотезу о линейности регрессии).

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

III.2.1. Формальная вероятностная основа контролируемого обучения

Для заданного распределения вероятностей g мы можем оценить плотность вероятности того, что функция f удовлетворяет уравнению (15). Эта плотность вероятности выражается как P( f |D, g), где D – вектор наблюдаемых данных обучения, и может быть вычислена с помощью формулы Байеса (18).

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

Математически это записывается как:

Аналогично,

Комбинируя уравнения (16) и (17) приходим к формуле Байеса:

В контексте обучения с учителем правило Байеса дает

Распределение вероятностей P( f |g) дает плотность вероятности того, что функция  f удовлетворяет уравнению (14) до получения информации, содержащейся в обучающем наборе данных. По этой причине оно известно как априорное распределение вероятностей. Распределение вероятностей P( f |D, g) представляет собой вероятность того же события после учета обучающих данных. Это известно как апостериорное распределение вероятностей. Распределение P(D| f, g), широко известное как функция правдоподобия, представляет собой вероятность наблюдения обучающих данных D при заданных f и g. Остаю- щийся член P(D|g) в формуле (18) не зависит от f и  может эффективно рассматриваться как нормировочная константа.

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

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

Подставляя (19) в (20), приходим к формуле:

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

Член в квадратных скобках в правой части уравнения (21) называется целевой функцией, которую нужно минимизировать. Функция, которая с максимальной апостериорной вероятностью удовлетворяет уравнению (15), приводит целевую функцию к минимуму. Идея поиска функции P(D| f, g), которая минимизирует целевую функцию, является общей для методов машинного обучения.

Первый член –lnP(D| f, g) в уравнении (21) представляет эмпирический риск – меру того, насколько хорошо данная функция воспроизводит обучающие данные. Предполагая, что все наблюдения в обучающих данных независимы, из уравнения (14) и определения g можно написать:

В более общем виде эмпирический риск записывают как:

III.2.2. Алгоритмы контролируемого обучения

Рассмотрим лишь несколько наиболее распространенных методов.

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

Тогда эмпирический риск принимает вид:

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

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

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

Выходным значением метода в этом случае является оптимальный набор w, определяемых в силу соотношения (27), как:

Уравнение (30) представляет собой тип аппроксимации методом наименьших квадратов, известный как регуляризация Тихонова [36]. Когда нормальная аппроксимация методом наименьших квадратов некорректна (например, когда величина XX^T сингулярна), регуляризация Тихонова может сделать задачу корректной.

б) Метод воспроизводящих ядер (kernel trick). Линейные методы наименьших квадратов очень популярны из-за высокой скорости и простоты, которые связаны с основным предположением, что выходные значения должно быть линейной функцией входных переменных. Во многих случаях это предположение нереально. Однако можно преобразовать исходное одномерное входное значение (x_i) в пространство признаков, обладающее большей размерностью (например, $({{x}_{i}},~x_{i}^{2},~x_{i}^{3}, \ldots ~x_{i}^{n})$), в котором уже возможна линейная регрессия [45, 46]. Их можно использовать в алгоритмах машинного обучения, просто подставив преобразованные входные переменные φ(x) вместо исходных входных переменных x. Однако многие алгоритмы, включая линейный метод наименьших квадратов, могут быть выражены способом, который зависит от входных переменных только через внутренние (скалярные) произведения, такие как (x₁ · x₂). В таких случаях достаточно найти функцию k(x_i, x_j), называемую ядром, которая позволяет определить скалярное произведение между φ(x_i) и φ(x_j):

Чтобы использовать алгоритм обучения в заданном пространстве функций, все, что нужно сделать, это использовать ядро интегрального преобразования k(x_i, x_j) вместо скалярного произведения (x_i · x_j) на протяжении всего алгоритма обучения. Такой способ в англоязычной литературе получил название как трюк с ядром (kernel trick) [46].

в) Гребневая регрессия c помощью ядер. Для воспроизводящего ядра k(x₁, x₂) рассмотрим пространство функций f, которые задаются формулой

Это функциональное пространство известно как гильбертово пространство воспроизводящего ядра [46].

В гильбертовом пространстве определяется норма ||   f   ||_H:

Воспроизводящие ядерные гильбертовы пространства обладают особыми свойствами, которые делают их особенно полезными для машинного обучения с учителем. Рассмотрим ситуацию, в которой пространство гипотез является гильбертовым пространством некоторого ядра. Если, например, используется гауссовское ядро, пространство гипотез будет состоять из линейных комбинаций гауссовых функций. Пусть член регуляризации задается как r(||   f   ||H), где r – монотонно возрастающая функция. Пусть эмпирический риск принимает общий вид $\mathop \sum \limits_i L\left( {{{y}_{i}},f({{x}_{i}})} \right)$, где сумма берется по всем точкам данных в обучающей выборке. Следовательно, целевая функция имеет вид:

Тогда теорема о представителе [47] утверждает, что функция, которая минимизирует эту целевую функцию, должна иметь вид:

Если используется квадратичная функция потерь и $r\left( {{\text{||}}\Delta \,f\,{\text{|}}{{{\text{|}}}_{{\text{H}}}}} \right) = {{\lambda ||}}\Delta \,f\,{\text{|}}{{{\text{|}}}_{{\text{H}}}}$ для некоторой положительной величины λ, то несложно показать, что целевая функция может быть записана как

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

г) Регрессия на основе гауссовского процесса (РГП) основана на несколько другом подходе в рамках байесовского описания [49]. Если в методах линейной и гребневой регрессии пытаются выразить выходную функцию f(x) через модельные функции путем подбора параметров разложения, то гауссовский процесс представляет f(x) косвенным, более тонким образом.

Гауссовским процессом называется случайный процесс $\xi (x,t)$, все конечномерные распределения которого для любого момента времени t нормальные, с ковариационной функцией C(x, x'). В случае машинного обучения гауссовский процесс – это набор случайных величин, любое конечное число которых имеет совместное гауссово распределение с заданным средним значением m и ковариационной функцией.

Рассмотрим задачу восстановления регрессии по обучающей выборке y_i = f(x_i), где величины y_i есть значения реализации случайного процесса (поля) при входных значениях x_i для i-го обучающего примера. Возникает задача прогноза значения поля у в новой точке x при условии, что в точках обучающей выборки поле наблюдались значения y_i:

Заметим, что по определению (40) для гауссовского случайного процесса c нулевым средним (m = 0)

Учитывая, что по определению условной вероятности

Метод РГП широко применяется в машинном обучении благодаря своему основному достоинству – отсутствию необходимости предугадывать вид функций, входящих в пространство гипотез. Вместо того чтобы утверждать, что f(x) относится к некоторым конкретным моделям (например, f(x) = mx + c), гауссовский процесс описывает эту функцию косвенным образом (непараметрическим), позволяя при этом обучающим данным четко “говорить” самим за себя. Однако, поскольку мы отказываемся определять выходную функцию f(x) относительно небольшим набором параметров, входящих в разложение по модельным функциям, то на каждом шагу регрессии приходится проводить вычисления по соотношению (42), куда входит весь набор обучающих данных. Это приводит к резкому росту вычислительных затрат при большом наборе тестовых примеров (кубический рост).

д) Искусственные нейронные сети прямого распространения. Использование нейронных сетей прямого распространения – это один из наиболее эффективных методов машинного обучения, активно развивающийся в последние годы [50]. Идея метода состоит в имитации работы биологического мозга посредством искусственных нейронных сетей (ИНС), о которых подробно сообщалось в предыдущем разделе. Пространство гипотез в нейронной сети определяется топологией соединений между узлами и параметризованными функциями активации, используемыми узлами. Оптимизация весов многослойной нейронной сети эффективно выполняется с помощью алгоритма обратного распространения, который эффективно минимизирует квадрат потерь с помощью градиентного спуска. Преимущества заключаются в высокой точности, высокой способности к обучению и высокой отказоустойчивости для зашумленных данных. Так же как и РГП, ИНС является непараметрическим методом, в котором заранее не делается предположений о виде функций f(x), в силу чего требуется относительно длительный процесс обучения. Это обстоятельство затрудняет использование метода “on the fly”, т.е. встраивание его в программы Монте-Карло или МД-моделирования.

III.3. Дескрипторы

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

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

Более широкую область применения и точность имеет метод кластерного разложения (cluster expansion, CE), который широко используется в различных приложениях для поиска материалов с требуемыми свойствами [53, 54]. Этот метод используется для расчета энергии конфигурации (но не сил, действующих на атомы) в структурах с кристаллической решеткой определенного типа. Поэтому часто потенциалы взаимодействия, построенные в модели СЕ, называют “потенциалами на решетке”. Кластерные разложения быстры и вполне успешны, когда стабильные структуры являются производными от определенной структуры (ГЦК, ОЦК и т.д.), но значительно менее полезны в более сложных случаях, когда атомные радиусы компонентов сильно отличаются и важными становятся эффекты релаксации структуры [55]. Мы не будем останавливаться на работах, использующих указанные подходы, поскольку они полно представлены в недавнем обзоре [56].

Далее сконцентрируемся на методах машинного обучения для регрессии потенциальной энергии как функции координат всех атомов системы, необходимых для проведения атомистического моделирования сложных соединений. В работах по моделированию конденсированных сред термин “дескриптор” употребляют реже, заменяя полный набор дескрипторов термином “структурное представление”. В данной работе мы не будем проводить различия между этими терминами, считая их эквивалентными. Мы знаем, что заряды ядер и положения атомов являются достаточными характеристиками для описания электронной плотности и полной энергии любой конденсированной системы [2], поскольку гамильтониан системы обычно полностью определяется этими величинами. В первых работах по конструированию MLIP для молекулярных систем [57] на вход нейронной сети подавались все декартовы координаты атомов системы. Вскоре были выявлены два основных недостатка такого подхода. Во-первых, большинство моделей машинного обучения требуют на входе фиксированного количества функций, поэтому переход от системы с одним количеством атомов в суперъ-ячейке к другой приводит к серьезным проблемам. Во-вторых, декартовы координаты плохо подходят в качестве дескрипторов, поскольку они не инвариантны относительно преобразований симметрии (вращения, трансляции, перестановки одинаковых атомов) суперъячейки сложного соединения, поэтому их число чрезвычайно избыточно. Действительно, суперъячейку одноэлементного кристалла можно задать всего тремя векторами трансляций.

В 2007 г. Бехлер и Парринелло [58] представили новый подход к решению вышеуказанных проблем (рис. 5). В их методе полная энергия E представлена в виде суммы вкладов от отдельных атомов E_i:

Рис. 5.

Структура нейросети, использованная в работе [58]. Декартовы координаты атома i обозначены радиус-вектором R_i. Координаты затем преобразуются в набор значений функций симметрии G_i, описывающих локальную геометрию окружения атома i, зависящую от расположения всех атомов в системе (дескриптор атомарного окружения). Значения дескриптора G_i подаются на вход нейронных подсетей S_i, на выходе которых формируется парциальный вклад E_i в полную энергию, связанный с i-м атомом.

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

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

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

Параметры ζ, η, λ и R_c являются подгоночными и определяются в процессе обучения нейронной сети.

Две указанные идеи, введенные впервые Бехлером и Паринелло [58], оказались настолько глубокими, что стали важным компонентом всех современных методов конструирования MLIP [59–65].

Действительно, представление энергии в виде суммы вкладов от отдельных узлов (43) очень близко по духу к идее потенциалов погруженного атома (EAM) о погружении нейтрального атома в “море” электронной плотности, создаваемой его ближайшим окружением [25]. При достаточно удачном выборе дескрипторов локального окружения можно добиться, что узельные энергии E_i, будут достаточно просто связаны с ними (по линейному, или близкому к нему закону). Тогда полная энергия E будет почти линейно зависеть от входных значений G_i , и мы приходим к простому случаю линейной регрессии.

Таким образом, любой метод конструирования MLIP по сути дела состоит из двух основных частей: 1) выбора дескриптора для описания локального окружения выбранного атома, 2) выбора эффективного алгоритма машинного обучения (регрессор), – реализующих регрессионную задачу в пространстве выбранных дескрипторов.

Нейронные сети были первой и самой популярной формой регрессора при использовании дескрипторов типа (45) и (46) для задачи конструирования MLIP [58–62]. В работах [63, 64] дескрипторы Бехлера и Паринелло уже использовались совместно с алгоритмом гребневой регрессией ядра.

Другой формой регрессора, предложенной в 2010 г., были гауссовские процессы, используемые в потенциалах гауссовой аппроксимации (GAP) [65]. При использовании с ядром гладкого перекрытия атомных позиций (SOAP) [66] они могут доказуемо аппроксимировать произвольное локальное многочастичное взаимодействие атомов, в отличие от потенциалов нейронной сети, используемых поверх двух- и трехчастичных дескрипторов [67]. Вероятно, третий по величине класс межатомных потенциалов основан на линейной регрессии с различными наборами дескрипторов локального окружения. Он включает в себя дескрипторы на основе разложения атомного кластера (Atomic cluster expansion) [68], и на основе спектрального анализа соседей (SNAP) [69, 70], включая недавнее расширение для многокомпонентных систем [71]. Одним из интенсивно развивающихся подходов к машинно-обучаемым межатомным потенциалам, появившимся в последние годы, является использование “потенциалов тензора моментов” (MTP), разработанных Шапеевым [72]. Расширение MTP на многокомпонентные системы [73] выходит за рамки линейной регрессии, однако существует альтернативная формулировка многокомпонентной MTP, которая остается линейной [74].

Существование столь большого числа дескрипторов локального окружения вызвало появление исследований [68, 75, 76], в которых проводится их сопоставление. В работе [77] было показано, что дескрипторы, которые изначально выглядят совершенно иначе, являются частными случаями общего подхода, в котором конечный набор базисных функций с возрастающими угловыми волновыми числами используется для разложения функции плотности в окрестности выбранного атома. В [68] было показано, что все перечисленные выше дескрипторы могут быть получены в рамках метода разложения атомного кластера, а в [75] была дана общая формулировка проблемы представления атомных структур в терминах (гладкой) плотности атомов.

Из вышеизложенного следует, что в настоящее время существуют 4 различных схемы машинного обучения межатомных потенциалов (MLIP), которые можно считать многообещающими. Это дескрипторы Бехлера и Паринелло (BP) в методе нейронных сетей, дескрипторы с ядром гладкого перекрытия атомных позиций (SOAP) вместе с регрессором, использующим гауссовские процессы с потенциалами гауссовой аппроксимации (GAP), метод на основе спектрального анализа соседей (SNAP) и метод “потенциалов тензора моментов” (MTP). Все 4 схемы показали очень хорошую точность, и в настоящее время не ясно, что лучше. Хотя схема BP кажется наиболее гибкой, GAP и MTP имеют лучшую производительность. В статье [72] было представлено сравнение между GAP и MTP и показано, что MTP быстрее и более точен для вольфрама. В недавней работе [78] было проведено всестороннее сравнение производительности GAP, MTP, BP и SNAP. С этой целью были использованы стандартизированные наборы DFT-данных для шести элементов (Li, Mo, Cu, Ni, Si и Ge), охватывающих различные кристаллические структуры (ГЦК, ОЦК и структура алмаза), химический состав и тип связи (металлические и ковалентные). Четыре MLIP оценивались с точки зрения их точности в воспроизведении энергий и сил DFT, а также свойств материала, таких как уравнения состояния, параметр решетки и упругие постоянные. Была предпринята попытка оценить требования к обучающему набору и относительные вычислительные затраты. В результате сравнения обнаружилось, что у всех использованных схем ключевые показатели (точность, объем обучающих данных и вычислительная стоимость) взаимосвязаны. Увеличение сложности схемы (т.е. количество степеней свободы), связанное с увеличением вычислительных затрат и увеличением обучающего набора, обычно приводит к более высокой точности. Для всех четырех схем построения MLIP существует “оптимальная” конфигурация, при которой дальнейшее увеличение количества степеней свободы дает лишь небольшое улучшение точности с увеличением вычислительных затрат. Было показано, что все схемы способны достигать точности, близкой к DFT, при прогнозировании энергии, сил и свойств материала, существенно превосходя традиционные IAP. Модели GAP и MTP демонстрируют наименьшие среднеквадратичные значения энергии и сил даже при весьма умеренном объеме обучающих данных. Однако модель GAP является одной из самых затратных с точки зрения объема вычислений для данной точности и продемонстрировала плохую экстраполяцию для полиморфных фаз в системах со структурой алмаза. Нейронносетевой подход BP показал худшую сходимость по энергии при небольших размерах обучающих данных, которая заметно улучшается при увеличении объема данных. Необходимо отметить несколько ограничений этой важной работы, которые являются возможными направлениями будущих исследований. Во-первых, сравнение проводилось для чистых элементов, не было предпринято никаких попыток включить в это исследование бинарные и, тем более, тройные системы. Во-вторых, остался открытым вопрос о том, что сильнее влияет на ключевые показатели – дескриптор или регрессор, поскольку не делалось попыток комбинировать различные дескрипторы локальной среды с различными структурами машинного обучения.

Обратим внимание читателей на новый дескриптор погруженной атомной плотности (EAD) в подходе нейросетевого машинного обучения [78], идея которого позаимствована из метода погруженного атома (EAM) [25]. В EAD электронная плотность заменяется квадратом линейной комбинации атомных орбитальных компонентов:

В заключение укажем на еще один недавно появившийся способ решения проблемы различной размерности входных данных при помощи графовых нейронных сетей [79].

IV.1. Высокоэнтропийные и мультиэлементные вещества

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

Ситуация резко изменилась после выхода в 2004 г. двух важных статей научных групп из Великобритании и Гонконга [80, 81], в которых предложили для разработки перспективных материалов изучать способы получения и свойства сплавов, состоящих из нескольких основных компонентов в равных или почти равных пропорциях. Мотивация статьи [80] состояла в необходимости изучить “неисследованную центральную область фазового пространства многокомпонентных концентрированных сплавов” (multiple principal elements alloys, MPEA). В работе [81] выдвигалась идея о том, что использование концентрированных смесей многих элементов в качестве основы нового сплава, приводящее к резкому увеличению конфигурационной энтропии, будет способствовать росту области однофазности системы при высоких температурах, блокируя выделения интерметаллидных фаз. В данном подходе сплавы с 5 или более основными элементами, концентрация которых составляет от 5 до 35 ат. %, получили название высокоэнтропийных сплавов, или материалов (high entropy alloys, HEA). Обе статьи сначала были встречены без особого внимания. Тем не менее к настоящему времени исследования в этой области увеличились до тысяч статей, поскольку во всем мире увидели огромный потенциал для поиска и разработки новых материалов с интересными и ценными новыми свойствами. Соответствующие ведомства США, Европейского Союза и Китая причислили HEA к ключевым направлениям развития на ближайшее будущее. HEA превосходят традиционные сплавы по механическим свойствам, коррозионной стойкости, стойкости к окислению и т.д., а часто и по стоимости, позволяя избегать использования дорогостоящих легирующих элементов. Продолжающийся поиск новых высокоэнтропийных систем (ВЭС) недавно расширился за пределы металлов и включил керамику со стабилизированной энтропией, такую как высокоэнтропийные оксиды [82], карбиды [83] и бориды [84].

В первые годы изучения HEA-материалов акцент в первую очередь делался на достижение высокой энтропии и, соответственно, на достижение максимального числа основных компонентов, что легче всего достичь в области металлических сплавов. Однако после подведения первых итогов масштабных исследований стало ясно, что исходная идея высокоэнтропийного сплава требует определенной корректировки [85]. В первых работах энтропия сплава оценивалась в приближении идеальных растворов, так что каждая из i компонент с концентрацией х_i создает вклад

Вопреки этому предположению термодинамический анализ большого числа НЕА показал, что идеальные растворы встречаются редко, поэтому простое увеличение числа основных компонентов не всегда приводит к существенному росту энтропии [86, 87]. Во-вторых, выяснилось, что рост числа основных элементов также существенно изменяет энтальпию. Причем из-за резкого увеличения количества различных пар соседствующих элементов возрастает вероятность того, что по крайней мере один тип таких пар приведет к изменению энтальпии, превосходящему рост энтропии.

Был сделан вывод, что формальное требование причислять к HEA лишь металлические системы, состоящие из 5 и более основных компонент, является слишком широким. Чтобы устранить этот барьер, был введен новый термин – комплексные концентрированные сплавы (ССА). Этот более всеобъемлющий термин сохраняет акцент на сложных по составу концентрированных сплавах, но включает в себя как HEA-сплавы, так и концентрированные тройные и четвертичные сплавы, допуская не только однофазные, но и многофазные растворы [88], включая высокоэнтропийные карбиды, нитриды и оксиды, а также объемно-аморфные материалы [89]. С учетом указанных различий в настоящем обзоре мы будем использовать термины «высокоэнтропийный материал» и «комплексный концентрированный материал» как синонимы.

IV.2. Прогнозирование стабильности высокоэнтропийных сплавов

Как уже отмечалось, основная идея НЕА состоит в том, чтобы за счет увеличения числа компонентов N максимизировать конфигурационную энтропию, чтобы добиться стабильности неупорядоченной фазы твердого раствора, подавляя таким образом образование интерметаллических фаз. Стабильность HEA-материала определяется термодинамическими правилами, вытекающими из основных принципов статистической физики Гиббса. Эти правила лежат в основе методологии CALPHAD [90], которая на протяжении многих десятилетий применяется для расчета фазовых диаграмм. Центральной величиной является экстенсивная энергия Гиббса всей системы $G\left( {P,T,\left\{ {{{n}_{i}}} \right\}} \right)$, являющаяся функцией давления P, температуры T и числа молей n_i для каждого элемента i = 1, …, N (N – количество элементов, для комплексных концентрированных систем N > 4). В силу экстенсивности свободная энергия может быть разложена на сумму существующих в системе фаз:

При проведении атомистического моделирования более удобно использовать молярные величины $\bar {G} = G{\text{/}}n$ и ${{\bar {G}}_{\alpha }} = {{G}_{\alpha }}{\text{/}}{{n}_{\alpha }}$, где n = $\mathop \sum \limits_i {{n}_{i}},$ n^α = $\mathop \sum \limits_i n_{i}^{\alpha }$, и мольные доли ${{x}_{i}} = {{n}_{i}}{\text{/}}n$ и $x_{i}^{\alpha } = n_{i}^{\alpha }{\text{/}}{{n}^{\alpha }}$. Тогда

Основная проблема феноменологических теорий – получить выражение для зависимости ${{G}_{\alpha }}(P,T,\{ x_{i}^{\alpha }\} )$. CALPHAD использует для этой цели базы данных термодинамической информации, основанные на экспериментальных измерениях (если таковые имеются), дополненные схемами численной интерполяции для заполнения тех областей фазовой диаграммы, где экспериментальная информация отсутствует.

Использование CALPHAD для изучения высокоэнтропийных систем сталкивается с двумя основными проблемами. Первая проблема чисто комбинаторная. Пусть проводится исследование системы из N = 6 компонент. Тогда только число 4-компонентных фазовых подобластей полного фазового пространства составит $\frac{{N!}}{{\left( {N~ - {\text{\;}}n} \right)!n!}} = 15$, что наглядно показывает объем требуемых экспериментальных данных. Вторая и основная проблема связана с доступностью и надежностью термодинамических баз данных [91]. Дело в том, что базы данных CALPHAD были разработаны для традиционных сплавов, поэтому полная термодинамическая оценка энергий Гиббса обычно дается только для составов, обогащенных одним основным элементом. Напротив, HEA и CCA лежат в центральных областях композиционного пространства, вдали от изученных границ, ограниченных бинарными и тройными подсистемами. Следовательно, вычисления в этих обширных пространствах композиций опираются на значительные экстраполяции, точность которых критически зависит от качества диаграмм систем более низкого порядка. Таким образом, надежность и точность прогнозов могут стать сомнительными, если они получены из неполных термодинамических описаний. Сенков и соавт. [92] определили критерии достоверности для расчетов CALPHAD, основанные на доле полностью термодинамически оцененных бинарных систем (F_AB) и доле полностью оцененных тройных систем (F_AT), включенных в базу данных. Они показали, что приемлемая точность прогноза достигается лишь в случае полной информации о всех бинарных системах и значительной доле тройных. Хотя база данных на основе Ni успешно предсказала равновесные фазы и фазовые доли в сплавах Al–Co–Cr–Cu–Fe–Ni [93], такие достижения остаются единичными. Действительно, даже специализированные базы данных для HEA, такие как TCHEA3 от ThermoCalc, содержат почти все описания бинарных систем для 26 элементов, но только 5% описаний тройных систем полностью оценены.

Отсюда следует, что предсказательные возможности метода CALPHAD существенно ограничены недостатком экспериментальных данных. Указанные проблемы вызвали появление полуэмпирических методов для предсказания того, какие многокомпонентные сплавы образуют твердый раствор. В них в качестве инструментов прогнозирования используются дескрипторы [94, 95], параметры которых соответствуют доступным, но ограниченным экспериментальным данным.

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

Совершенно ясно, что огромное композиционное пространство, охватываемое HEA и CCA, не может быть исследовано только экспериментальными усилиями. Необходимо подключить все возможности компьютерного моделирования. Особенно важным вычислительным инструментом являются вычисления, основанные на первых принципах (аb initio). Такие расчеты проводятся в рамках теории функционала плотности (DFT), основаны только на квантово-механических законах и физических константах и, таким образом, позволяют предсказывать и исследовать свойства материалов без эмпирических данных. Далее мы позволим себе использовать термины ab initio и DFT как взаимозаменяемые.

Одним из магистральных путей такого моделирования является объединение методов CALPHAD с ab initio расчетами полной энергии вещества, чтобы использовать относительные преимущества каждого из них. Первопринципные расчеты превосходно подходят для предсказания низкотемпературных значений энтальпий, которые часто просто недоступны экспериментально, а также параметров структуры и магнитных моментов. CALPHAD отлично подходит для моделирования температур переходов и топологии фазовой диаграммы. Теория функционала плотности [1, 2], на которой основано ab initio моделирование, в своей первоначальной формулировке является теорией основного состояния, т.е. строго позволяет рассчитать только энергию данной конфигурации атомов в суперъячейке при 0 К. Однако энтропийный эффект, присущий высокоэнтропийным материалам, проявляется только при высоких температурах, где важную роль играют энергии возбуждений над основным состоянием – фононов, магнонов и т.д. Поэтому необходимо комбинировать DFT-вычисления с термодинамическими концепциями и статистическим усреднением.

Такая комбинация позволяет рассчитать свободную энергию Гельмгольца ${{F}_{\alpha }}\left( {V,T,\left\{ {{{n}_{i}}} \right\}} \right)$ при конечных температурах без использования данных эксперимента [96]. После этого величину G_α получают с помощью стандартного соотношения:

Выражение (53) показывает, что свободная энергия зависит не только от типа решетки фазы α, но и от типа ближнего порядка в ней. Учет ближнего порядка при проведении DFT-расчетов достаточно сложен, поэтому мы не будем его рассматривать, отправляя читателей к обзорам [96, 98].

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

Здесь ${{E}_{\alpha }}(V,\{ x_{i}^{\alpha }\} )$ – полная электронная энергия атомной структуры α-фазы, параметр решетки которой и положение всех атомов отрелаксированы с помощью DFT для достижения минимума при T = 0 K; $S_{\alpha }^{{{\text{conf}}}}(\{ x_{i}^{\alpha }\} )$ – конфигурационная энтропия, $F_{\alpha }^{{{\text{el}}}}(V,T,\{ x_{i}^{\alpha }\} )$ – свободная энергия, обусловленная электронными возбуждениями, $F_{\alpha }^{{{\text{mag}}}}(V,T,\{ x_{i}^{\alpha }\} )$ – свободная магнитная энергия, а $F_{\alpha }^{{{\text{vib}}}}(V,T,\{ x_{i}^{\alpha }\} )$ – свободная энергия, обусловленная колебаниями атомов.

Все указанные вклады вполне поддаются DFT-расчету [96]. Неупорядоченная химическая конфигурация моделируется на основе большой, но конечной суперъячейки с использованием периодических граничных условий. Большая суперъ-ячейка приводит к значительным вычислительным затратам, поэтому на практике расчеты для HEA-материалов выполняются с использованием базиса плоских волн и метода проекции присоединенной волны (PAW), реализованного в пакете моделирования VASP [99, 100]. В принципе, любая атомная конфигурация, полученная путем однородного случайного распределения, может быть использована для моделирования химического беспорядка внутри суперъячейки. Однако, чтобы достаточно хорошо представить идеальную химическую хаотичность, соответствующие суперъячейки потребуют нескольких сотен атомов. Поэтому более удобным является подход SQS [101], в котором частицы в суперъячейке умеренного размера (несколько десятков узлов) распределяются таким образом, чтобы минимизировать функции корреляции пар атомов. Набор полученных таким образом результатов вполне достаточен, чтобы дальнейшее уточнение ближнего порядка в системе проводить методами решеточного Монте-Карло [102].

Заметим, что моделирование HEA-материалов из первых принципов сильно осложняется большим количеством компонент, что влечет за собой комбинаторно высокие вычислительные затраты. Дополнительное усложнение в ряде случаев привносит необходимость учета беспорядка в сложной химической и магнитной структуре [94, 96]. Поэтому характерная особенность ab initio моделирования HEA-материалов состоит в обязательном использовании при вычислении обменно-корреляционного функционала электронной плотности приближения обобщенного градиента(GGA). Это связано с тем, что большинство интересных HEA-систем являются магнитными и, следовательно, требуют использования спин-поляризованной версии DFT [103]. Как известно, основным успехом GGA было правильное предсказание ферромагнитного ОЦК-состояния как основного для железа [104]. По этой причине для точного определения фазовой стабильности систем, содержащих магнитоактивные компоненты Fe, Cr, Ni, использование приближения GGA является необходимым. К настоящему времени проведены DFT-расчеты более 2000 HEA-материалов, информация о которых приведена в недавних обзорах [96, 105, 106].

IV.3. Методика атомистического моделирования HEA с применением MLIP

Первоначальное изучение фазовых равновесий в комплексных концентрированных соединениях проводилось методами вычислительной термодинамики (CALPHAD), которые позволяют экстраполировать фазовые диаграммы высокого порядка из термодинамических описаний систем низшего порядка (т.е. от унарных к двойным, от двойных к тройным и т.д., вплоть до n-го порядка). Рассмотренные в предыдущей главе ограничения существенно снижают потенциальную точность прогнозов CALPHAD применительно к высокоэнтропийным системам. С другой стороны, методы ab initio моделирования, о которых говорилось выше, имеют принципиальные ограничения по пространственно-временному масштабу моделирования [96]. Поэтому возможность получать с помощью методов машинного обучения потенциалов межатомных взаимодействий с точностью, близкой к точности DFT-расчетов, в случае высокоэнтропийных систем является прорывной идеей. Действительно, использование таких потенциалов в методах Монте-Карло и молекулярной динамики практически снимает все ограничения на моделирование механических, диффузионных и электронных характеристик.

Раннее атомистическое моделирование HEA проводилось либо с полуэмпирическими EAM-потенциалами, либо базировалось на совместном использовании DFT-расчетов в приближении когерентного потенциала (CPA) [107]. Данный подход в сочетании с моделированием Монте-Карло (MC) широко применялся в последние несколько лет для решения проблемы фазовой стабильности и ближнего порядка многокомпонентных сплавов [108–110] из-за высокой вычислительной эффективности. Однако использование метода среднего поля, характерное для когерентного потенциала, не в состоянии учитывать локальные релаксационные эффекты, которые могут быть важны для рассмотрения фазовой стабильности [110]. Более точным является моделирование методом Монте-Карло на основе DFT с использованием суперъячейки [111]. Эти вычисления, однако, чрезвычайно требовательны в вычислительном отношении и могут применяться лишь для сплавов с числом компонент не более четырех. Поэтому чаще всего данные о полной энергии ограниченного набора структур, полученных с помощью DFT, отображают на эффективных моделях, таких как метод кластерного разложения [112]. Однако последние данные показывают (см., например, обзор [113]), что использование кластерного разложения, равно как AIMD, является вычислительно неэффективным для многокомпонентных сплавов. Подробные обзоры атомистического моделирования HEA такого рода представлены в обзорах [96, 114].

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

Нейронные сети были первым методом машинного обучения при построении поверхностей потенциальной энергии. Одним из наиболее успешных приложений машинного обучения для создания надежного представления поверхности потенциальной энергии является подход Белера и Парринело [58], в котором впервые было предложено представлять полную энергию системы как сумму атомных энергий:

С момента публикации в 2007 г. в подход Белера и Парринело (NNP) было внесено несколько улучшений. В статье [115] был предложен метод учета дальнодействующих кулоновских взаимодействий для ионных систем путем уравновешивания зарядов через нейронные сети, что позволило воспроизвести несколько объемных свойств CaF₂. Затем в статье [116] разработана стратегия обучения иерархических многокомпонентных систем, последовательно переходя от элементов к бинарным, а затем к тройным и т.д. системам. Тем не менее метод NNP оказался не слишком удобен для многокомпонентных систем с числом компонент N > 3. Это связано с использованием в качестве дескриптора атомно-центрированных функций симметрии. Для многокомпонентных систем функции симметрии должны быть определены отдельно для разных пар (радиальная часть) и троек атомов (угловая зависимость). В результате сложность метода экспоненциально возрастает с ростом N [117, 118]. Хотя существуют примеры использования NNP для четырехэлементных [119] и пятиэлементных [78] систем, их точность, как указано в [117], значительно ниже, чем у NNP, построенных для систем с одним-тремя элементами.

Более перспективными для многокомпонентных систем оказались методы гауссовских процессов с потенциалами гауссовской аппроксимации (GAP) [65] при использовании с ядром гладкого перекрытия атомных позиций (SOAP) [66]. Во-первых, они позволяют аппроксимировать произвольное локальное многочастичное взаимодействие атомов, в отличие от потенциалов нейронной сети, использующих лишь двух- и трехчастичные дескрипторы. Во-вторых, в методе GAP для многокомпонентных систем можно использовать весовые коэффициенты для каждого типа атомов при построении плотности соседей, поэтому вычислительные затраты растут значительно медленнее, чем у NNP [120]. Метод SNAP [69, 70] также был обобщен на многокомпонентные системы [71] и c успехом использовался для моделирования системы Be–W [121].

Вполне пригодным оказался и разработанный Шапеевым подход потенциалов тензора моментов (МТР). Метод легко различает вклад различных пар частиц в радиальные функции [73], поэтому количество параметров модели растет с увеличением числа компонент менее чем квадратично. В недавней работе [122] было проведено объективное сравнение подходов GAP-SOAP и МТР при моделировании с помощью MLIP-энергий конфигураций сплава Ag–Pd в широком диапазоне составов. Оба типа потенциалов обеспечивают превосходную точность описания фононных спектров, близкую к точности DFT, а также точек плавления и фазовых переходов. При этом метод МТР показал очень хороший баланс между точностью и вычислительной эффективностью.

В работе [89] Драутц недавно развил идею кластерного разложения применительно к построению дескрипторов и получил новый класс MLP, названный “разложением по атомным кластерам”. Основываясь на этой идее, Шапеев разработал малоранговый потенциал межатомного взаимодействия на решетке (LRP) [123]. Как и другие потенциалы межатомного взаимодействия, LRP использует соотношение (10), т.е. полная энергия системы представляется в виде суммы вкладов отдельных атомных окружений E_i. Далее предполагается, что атомы исследуемой системы не сильно смещены относительно узлов идеальной решетки, так что координационное число N_c ближайших соседей в пределах радиуса обрезания R_с не изменяется. Тогда вклад i-го атома в энергию системы однозначно определяется типом атома i и набором типов атомов (Fe, Ti, Mn и т.д.), входящих в его ближайшее окружение. Следовательно, E_i = V(σ₁, . . . σ_n), где V(σ₁, . . . σ_n) – потенциал межатомного взаимодействия, σ₁, . . . σ_n – набор, характеризующий тип атома i и типы его ближайших соседей (n = Nc + 1). Потенциал V(σ₁, . . . σ_n) с достаточной точностью представим в виде произведения малоранговых тензоров:

IV.4. Активное обучение

Известно, что MLIP дают плохое описание (большие ошибки подгонки) для областей PES, где обучающих точек недостаточно, как показано на рис. 6. Из-за огромного числа локальных минимумов на PES невозможно построить обучающий набор, охватывающий все локальные области, поэтому никогда нет гарантии, что MLIP предоставляет разумные описания для всех низкоэнергетических потенциальных ям. В наиболее серьезных случаях, когда MLIP дает плохое описание областей PES, связанных с глобальным минимумом, поиск структуры завершится ошибкой. Одним из решений для преодоления этой проблемы являются схемы активного обучения, использующие идею обучения на лету (on the fly) [126, 127]. Суть идеи – доверить уточнение обучающей выборки при проведении итераций компьютеру. Компьютерная программа определяет, находится ли возникшая при итерации структура за пределами эталонного набора данных. Если это так, то для нее проводится разовый DFT-расчет и полученные данные добавляются к обучающему набору. Таким образом, MLIP обновляется итеративно во время поиска структуры путем постоянного включения экстраполяционных конфигураций в обучающий набор.

Рис. 6.

Схематическая иллюстрация поверхности потенциальной энергии соединения в многомерном пространстве конфигураций, а также связь между референсными данными DFT-расчетов и MLIP. Сплошная линия обозначает реальный ход PES, а пунктирная линия обозначает подогнанный PES с помощью MLP. Серым фоном выделена область с недостаточным числом точек тренировочного набора. Как правило, MLP отлично подходят для регионов, плотно охваченных контрольными тренировочными точками (черные точки), и наоборот.

Таким образом, в настоящее время для моделирования комплексных концентрированных соединений на основе MLIP-потенциалов используется итерационный многомасштабный подход (рис. 7), который включает несколько основных этапов. На первом этапе выбирается решетка моделируемого вещества, узлы которой случайным образом заполняются компонентами (от 10 до 90 ат. % для каждого компонента) и таким образом генерируются начальные приближения для нескольких сотен репрезентативных структур. На втором этапе производится DFT-моделирование указанных репрезентативных структур, в результате которого определяются равновесные значения параметра решетки и полной энергии системы. Полученные результаты объединяются в базу данных. На следующем шаге на основе базы квантово-механических расчетов производится параметризация MLIP-потенциала на основе одного из методов, рассмотренных нами выше в разделе II. Далее полученный потенциал используется для моделирования эволюции атомной структуры методами Монте-Карло или молекулярной динамики (МД). По результатам моделирования отбираются наиболее репрезентативные конфигурации ВЭС. Затем для них определяется ошибка представления полной энергии по отношению к результатам DFT-расчетов. Если ошибка превосходит 5 мэВ/атом, то цикл самосогласования, изображенный на рис. 7, переходит на следующий виток. А именно, база данных дополняется новым набором репрезентативных конфигураций, на обновленной базе происходит обучение потенциала и т.д.

Рис. 7.

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

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

Алгоритм активного обучения на лету был адаптирован для оптимизации системы потенциалов нейронной сети (NNP) Артитом и Белером [128] и включен в пакет ÆNET [129], реализующий метод Белера–Паринелло. В последние годы схема итеративного построения MLIP была предложена Дерингером для метода гауссовской аппроксимации GAP [130] и Шапеевым и соавт. для потенциалов тензора моментов МТР [131]. Все эти приложения показали, что обучение и производственный цикл можно ускорить на несколько порядков.

V.1 Фазовые диаграммы

Область, в которой использование потенциалов машинного обучения принесло значительный прогресс, – это расчет термодинамических свойств, таких как фазовая стабильность и фазовые диаграммы, а также другие величины, зависящие от температуры. Основная причина успеха MLIP для термодинамических предсказаний HEA – это локализация взаимодействий в системе при в высоких температурах: дальнодействующие при 0 K межатомные взаимодействия, имеющие электронную природу, существенно ослабляются, когда возникают атомные колебания и нарушается симметрия кристалла [132]. В этом случае вычисление свободной энергии при высокой температуре представляет задачу, которая идеально подходит для потенциалов, характеризующихся ограниченным радиусом действия. Превосходная точность методов машинного обучения потенциалов для HEA была продемонстрирована в недавней работе [133], где МТР-метод Шапеева был использован для вычисления свободной энергии теплостойких ОЦК-сплавов системы NbMoTaVW при 3000 K с точностью, близкой к DFT-результатам. Численно точные колебательные свободные энергии были получены двухэтапным термодинамическим интегрированием с использованием метода динамики Ланжевена [134].

Ряд других примеров также подтверждает успешные возможности машинного обучения для анализа высокотемпературных свойств HEA: метод потенциалов низкого ранга (low range potentials, LRP) [122] совместно с Монте-Карло-моделированием был применен для исследования фазовых диаграмм ОЦК MoNbTaW и ГЦК CrFeCoNi [124, 135], причем в обоих случаях были обнаружены ранее неизвестные структуры. В работе [136] потенциал нейронной сети был использован для исследования структуры тугоплавкого высокоэнтропийного расплава.

Возможности MLIP-потенциалов позволяют рассчитывать полную свободную энергию системы при различных температурах и объемах, включая точные ангармонические составляющие. Исходя из этого могут быть точно рассчитаны такие свойства материала, как объемный модуль, теплоемкость и тепловое расширение. Хорошим примером может служить работа [137], в которой рассмотрена среднеэнтропийная тройная ГЦК-система FeCoNi и были рассчитаны локальные искажения решетки и упругие постоянные для ряда температур.

V.2. Моделирование влияния ближнего порядка на свойства HEA

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

Ближний порядок обычно количественно описывается параметрами Уоррена–Каули [138], которые определяются парными корреляциями корреляций между элементами i и j:

БП является локальной структурной характеристикой, которую трудно зафиксировать экспериментально. Чтобы компенсировать этот недостаток, для ее расчета в HEA были успешно применены различные атомистические методы. Наиболее популярные подходы включают моделирование методом Монте-Карло с использованием потенциалов машинного обучения [124]. В качестве примера на рис. 8 показана температурная зависимость параметров Уоррена–Каули для первой оболочки четырех различных пар в четырехэлементном HEA сплаве MoNbTaW [125].

Рис. 8.

Зависимость параметров ближнего порядка Уоррена–Каули от температуры в MoNbTaW для четырех различных пар атомов по данным Костюченко [125].

Видно, что для этого сплава две фазы B2 (Mo, Ta) и B32 (Nb, W) разделяются при низкой температуре, в то время как фаза B2 (Mo, W; Nb, Ta) с дальним порядком стабильна примерно до 600 К. Следует отметить также, что параметры Уоррена–Каули не равны нулю даже при относительно высокой температуре, что свидетельствует о наличии БП даже для неупорядоченного твердого раствора.

Ближний порядок определяется не только особенностями химического взаимодействия компонентов. Для сплавов, содержащих 3d-элементы, БП во многом определяется их магнитным взаимодействием. Так, зависимость между БП и магнетизмом имеет решающее значение для сплавов семейства CrMnFeCoNi. Это семейство представляет очень своеобразную магнитную ситуацию, в которой ферромагнитные металлы (Fe, Co и Ni) и антиферромагнитные металлы (Cr и Mn) перемешаны случайным образом. Благодаря этому конкурирующие магнитные обменные взаимодействия между атомными спинами с разной ориентацией могут приводить к эффектам магнитной фрустрации, когда, например, атомы хрома вынуждены отталкиваться друг от друга, поскольку они имеют параллельные спины из-за ферромагнитного окружения. Приведенное соображение стало основой для нескольких исследований БП в тройных и четверных сплавах CrCoNi и CrCoFeNi соответственно [139, 140]. И экспериментальное, и теоретическое исследования показали, что параметры Уоррена–Коули существенно отличаются от случайного твердого раствора: среднее количество пар ближайших соседей Cr–Cr было уменьшено примерно на 40%, тогда как количество пар Ni–Cr было увеличено. Необходимая энергетическая движущая сила для БП в этом сплаве обеспечивается магнитными обменными взаимодействиями.

Следовательно, в сплавах CrCoFeNi магнитные моменты Fe, Co и Ni в основном ориентируются параллельно, в то время как атомы Cr предпочитают антиферромагнитное упорядочение по отношению к ним.

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

V.3. Ближний порядок и механические свойства ВЭС

Несмотря на свой локальный характер, ближний порядок может существенно влиять на макроскопические свойства ВЭС. В последние годы эксперименты позволили установить, что БП глубоко влияет также на механические свойства ВЭС. Так, в работе [141] было проведено ab initio моделирование высокоэнтропийных сплавов систем TiZrHfTa и TiZrNbHf. Данные ОЦК-сплавы демонстрируют превосходные свойства, такие как высокая прочность и высокая температура плавления, и, таким образом, являются потенциальной альтернативой суперсплавам на основе никеля. Уникальной механической прочности данных ВЭС, связанной с локальными искажениями структуры, однако, сопутствует низкая пластичность, что характерно для сплавов с ОЦК-решеткой. Как показал эксперимент, можно подобрать состав системы, чтобы достичь высокой пластичности, индуцированной трансформацией (TRIP), при фазовом переходе между ОЦК- и ω-фазами. Исследование показало, что ближний порядок и связанные с ним локальные искажения решетки в значительной степени влияют на фазовую стабильность изученных ВЭС. При правильном учете релаксационных эффектов предсказанные составы вблизи энергетического равновесия ОЦК–ГПУ оказались близкими к экспериментальным составам, для которых наблюдаются хорошая прочность и пластичность за счет TRIP-эффекта.

Ближний порядок может кардинально изменить механические свойства ГЦК-HEA. Было обнаружено, что различные локальные химические структуры влияют на энергию дефекта упаковки (ЭДУ) в тройных, четверных и пятикомпонентных сплавах семейства CrMnFeCoNi [142, 143] до такой степени, что ЭДУ даже качественно изменяется от отрицательного к положительному значению. Более высокие степени БП обычно увеличивают ЭДУ, потому что для создания дефекта упаковки необходимо локально нарушить порядок в системе, хотя бывают и исключения. Увеличение ЭДУ в ГЦК-ВЭС может изменить доминирующее деформационное поведение от мартенситного превращения к образованию двойников или простому скольжению дислокаций.

V.4 Новые типы дефектов в ВЭС: дефекты внедрения и поверхности

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

Влияние легирования углеродом в ВЭС семейства материалов с элементами Cr, Mn, Fe, Co и Ni изучалось в экспериментальной работе [144], где было показано, что добавление этого элемента увеличивало текучесть и предел прочности, а также влияло на пластичность сплава. В дальнейшем эта же система исследовалась с помощью ab initio моделирования [145]. Анализировалось влияние межузельных атомов углерода на CrMnFeCoNi в ГЦК- и ГПУ-фазах. Было обнаружено, что добавление углерода энергетически стабилизирует ГЦК-фазу по отношению к ГПУ-фазе и, таким образом, увеличивает ЭДУ в ГЦК-фазе системы CrMnFeCoNi, что согласуется с экспериментальным наблюдением [146]. Еще одно применение межузельного легирования в ВЭС относится к хранению водорода. Hu и соавт. [147] исследовали TiZrNbMoHf с помощью расчетов из первых принципов и обнаружили фазовый переход ОЦК–ГЦК при высоком содержании атомов водорода, что согласуется с экспериментальными данными [148]. Они также обсудили связь между растворимостью водорода и искажениями решетки и сообщили, что искажения решетки в ВЭС может повысить растворимость водорода. На данном этапе атомистическое моделирование ВЭС как потенциальных кандидатов для хранения водорода находится на начальном этапе, однако можно надеяться на интересные результаты.

V.5. Поверхность вещества

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

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

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

Комбинация расчетов из первых принципов и машинного обучения для пятикомпонентных HEA IrPdPtRhRu, CoCuGaNiZn и AgAuCuPdPt [149, 150] была использована для поиска новых катализаторов для реакций восстановления O₂, CO₂ и CO: для этих сплавов база данных энергий поглощения для различных перестановок элементов, окружающих поглощенную молекулу, была приспособлена с использованием в качестве дескрипторов химических признаков атомов в первых трех соседних оболочках вокруг абсорбента. Полученная модель затем использовалась для определения наилучшего состава путем максимального увеличения количества участков с заранее определенной оптимальной энергией поглощения, что привело к значительному увеличению расчетной каталитической активности. По мнению авторов ВЭС-катализаторы выглядят весьма перспективным направлением исследований.

V.6. Высокоэнтропийная керамика

После открытия в 2015 г. высокоэнтропийных оксидов, в которых однозначно был продемонстрирован энтропийный переход к однофазной системе [81], круг высокоэнтропийных материалов стал быстро расширяться. К ним стали причислять бориды [83], карбиды, нитриды, сульфиды и силициды высокоэнтропийных сплавов, которые обладают повышенными свойствами, ведущими к их практическому применению в теплоизоляции и защите от коррозии, термоэлектричестве, расщеплении воды, катализе и накоплении энергии [105]. Сложные механизмы, лежащие в основе этих свойств, которые часто трудно установить, становятся перспективными возможностями для вычислительного моделирования. Важным шагом в этом отношении явилась работа [151], в которой стандартный DFT-подход был существенно видоизменен для расчета энергии образования многокомпонентных систем с различными типами химической связи, что существенно облегчило количественный анализ. Так, в работе [152] с помощью теории функционала плотности в сочетании со специальными квазислучайными структурами (SQS) были изучены структурные, механические и электронные свойства высокоэнтропийного карбида (TaNbHfTiZr)C в диапазоне давлений 0–50 ГПа. Было показано, что пластичность (TaNbHfTiZr)C значительно улучшается с увеличением давления, а переход от хрупкого характера разрушения к пластичному происходит примерно при 20 ГПа.