Молекулярно-динамическое моделирование стабильности гидридов титана при помощи машинно-обученных потенциалов

Титан и его сплавы рассматриваются как перспективные материалы для хранения водорода [1]. Кроме того, определенную актуальность изучение гидридов титана имеет с точки зрения водородного охрупчивания конструкционных материалов. Согласно экспериментальным данным [2], со стехиометрическим составом стабильными являются гидриды с композицией TiH 2 . Фаза δ -TiH 2 имеет кристаллическую структуру типа флюорита с пространственной группой Fm3m , в которой атомы титана располагаются в гранецентрированной кубической решетке (ГЦК), а водород занимает тетраэдриэческие междоузлия решетки. При температурах ниже ~310 K происходит [3, 4] превращение в гранецентрированную тетрагональную ε -фазу. Кубическая δ -TiH x фаза также наблюдается и для нестихиометрического содержания водорода при 1,5 <  x <  2,0. Причем превращение в тетрагональную фазу при комнатной температуре 298 K, согласно эксперименту [4], происходит при содержании водорода выше x = 1,84. Кроме того, дискуссионным является также вопрос, является ли фаза γ -TiH c соотношением c/a >  1 метастабильной или термодинамически устойчивой. Данная структура была обнаружена в некоторых гидрированных образцах (от 1 до 3 ат. % водорода) α -титана (ГПУ-решетка) в виде выделений в α -матрице[5]. Первопринципные расчеты свойств гидридов титана выполнялось неоднократно. Изучались как профили энергий относительно тетрагональной деформации при нулевой температуре [6], так и термодинамические свойства в рамках фононного подхода [7, 8]. В работе [9] в рамках квазигармонического приближения температура ГЦК-ГЦТ превращения для TiH 2 была определена выше 800 К. Учет электронного вклада в свободную энергию Гиббса при конечных температурах, проведенный в работе [10], позволил получить ε - δ переход при 418 K. Фазовые превращения между кубической и тетрагональной фазами с использованием метода молекулярной динамики (МД) для мартенсита Fe-C исследовалось ранее в работах [11–13].

Физика

В последнее десятилетие бурное развитие получил подход по разработке полуэмпирических потенциалов для молекулярно-динамического моделирования с использованием алгоритмов машинного обучения и машинно-обученных межатомных потенциалов (MLIP). В рамках такого подхода потенциалы для описания межатомных взаимодействий строятся не в виде физически обоснованной функциональной зависимости со специально подобранными параметрами, а в виде набора (часто весьма большого) простых с вычислительной точки зрения дескрипторов. Построенные таким образом потенциалы с помощью алгоритмов машинного обучения параметризуются с использованием результатов первопринципных расчетов энергетических и силовых характеристик моделируемых ансамблей и отдельных частиц. В работе [14] было проведено построение машинно-обученного потенциала для системы Ti-H и проведен расчет диффузии и механических свойств гидридов титана с различным содержанием водорода. Работа [15], в которой было выполнено обучение межатомного потенциала с помощью нейросетей, посвящена изучению формирования гидридов с ГЦК решеткой и её кристаллографической связью с матрицей α -Ti. Однако, в обеих работах методике обучения потенциала уделено недостаточно внимания, например, нет сведений ни о размере обучающей выборки, ни о радиусе обрезания потенциала.

С учетом сказанного выше, представляется актуальным построение MLIP-потенциалов для гидридов Ti-H с ГЦК решеткой и МД-моделирование их поведения при конечных температурах для наблюдения ГЦК-ГЦТ превращений и влияния на них содержания водорода.

Обучение машинно-обученного потенциала

Ab initio моделирование титана и его гидридов в данной работе выполнялось в рамках теории функционала плотности (DFT) с псевдопотенциальным подходом. Использовался метод Projector Augmented Wave (PAW) [16], реализованный в пакете Vienna Ab Initio Simulation Package (VASP) [17-19]. Обменно-корреляционное взаимодействие рассматривалось с использованием обобщенного градиентного приближения (GGA) в формулировке Perdew–Burke–Ernzerhof (PBE) [20]. В рассматриваемом псевдопотенциале для титана учитывалось 12 внешних электронов в качестве валентных. Размытие электронных состояний было выбрано в форме Methfessel–Paxton [21] c шириной σ = 0,05 эВ. Энергии обрезания плоских волн задавалась на уровне 600 эВ. Разбиение первой зоны Бриллюэна в обратном пространстве было задано k -сеткой 4×4×4 с центром в точке Γ, что обеспечивает максимальное расстояние между точками на уровне 0,15 Å^–1.

В данной работе в качестве MLIP модели использовался Moment Tensor Potential (MTP) [22– 24] реализованный в пакете MLIP-3 [25] и интегрированный в пакет классической молекулярной динамики LAMMPS [26]. Процедура обучения потенциала заключалась в генерации некоторого количества снимков, взятых из первопринципных молекулярно-динамических траекторий с интервалом в 100 шагов по 1 фс, чтобы избежать корреляций между выбранными конфигурациями. Для этого выполнялись МД расчеты в пакете VASP для систем TiH x ( x = 2,0, 1,75, 1,5, 1,0) с различными степенями тетрагональности c/a от 0,95 до 1,05. Размер системы составлял 108 атомов титана, расположенных в узлах ГЦК решетки, а атомы водорода занимали тетраэдрические междоузлия. Причем системы с нестехиометрическим соотношением создавались путем удаления случайных атомов H в решетке дигидрида титана. Также для создания обучающей выборки выполнялись МД расчеты с различными соотношениями c/a для упорядоченной фазы γ -TiH, а также для системы TiH, в которых водород находился в октаэдричских порах. Всего таким образом было получено 1020 снимков.

Затем мы использовали процедуру активного обучения (Active Learning, AL), подробно представленную в работе [27], которая позволяет эффективным образом покрывать конфигурационное пространство дескрипторов потенциала и, таким образом, проводить дообучение модели только на необходимых конфигурациях. После обучения MTP на основе начальных данных выполнялись МД расчеты, в ходе которых на каждом временном шаге происходило вычисление степени экстраполяции γ на основе подхода D- оптимальности для получаемых атомных конфигураций, которая показывает есть ли в текущей обучающей выборке конфигурация «похожая» на текущее расположение атомов. Если степень экстраполяции γ находилась в области 2 <  γ < 10, то, такие конфигурации добавлялись к обучающему набору. Если же γ становилась больше, чем γ cr = 10,0, то происходила остановка расчета, поскольку в этом случае высока вероятность попасть в нефизичные области конфигурационного пространства. Эти выбранные снимки исполь-

Чирков П.В., Молекулярно-динамическое моделирование стабильности Мирзоев А.А. гидридов титана при помощи машинно-обученных потенциалов зовались в DFT для расчета сил, энергии и напряжений. Результаты добавляли в обучающие наборы MTP, после чего происходила дооптимизация коэффициентов MTP, далее МД расчеты продолжались. Итеративная процедура активного обучения выполнялась до тех пор, пока для систем TiH x ( x = 2,0, 1,75, 1,5, 1,0) на протяжении 100 тыс МД шагов в расчете не перестанут появляться конфигурации с γ > 2,0. По мере добавления новых конфигураций в обучающую выборку и улучшения качества работы потенциала температура в расчетах увеличивалась от 300 до 1000 К.

В итоге общее количество снимков при обучении составило 2963 конфигурации и, соответственно, такое же количество значений потенциальных энергий. Количество рассчитанных компонент сил, действующих на атомы титана составило 3·2963·108 ~ 9 ∙ 10⁵, количество компонент сил, действующих на атомы водорода, менялось в зависимости от концентрации, но имело примерно такой же порядок. Кроме того, для оптимизации потенциала использовались значения тензора напряжений (только его вириальная часть). Причем относительные веса для энергий, сил и давлений равнялись 1, 0,25, 0,25 соответственно. Радиус обрезания потенциала был выбран на уровне 5,0 Å, причем для больших значений r cut вплоть до 6,0 Å не было отмечено увеличение точности оптимизации MTP, что свидетельствует о короткодействующем характере взаимодействия в гидридах титана. Радиус обрезания r cut = 5,0 Å обеспечивает рассмотрение для ГЦТ-фазы TiH 2 74 ближайших соседей, из которых 18 атомов являются титаном. Для изучения свойств гидридов титана с различным содержанием водорода были обучены MTP различных уровней, для которых в таблице приведены значения ошибок для энергий и сил на обучающем наборе данных. Для дальнейших расчетов нами был выбран потенциал 22 уровня, содержащий 421 числовых коэффициента, определенных в ходе оптимизации.

Среднеквадратичные ошибки (СКО) по силам ( F ) и энергии ( E ) при обучении MTP потенциалов различного уровня

«уровень» MTP	кол-во коэффициентов	СКО E , мэВ/атом	СКО F , мэВ/Å
16	130	13,6	163
18	201	8,42	127
20	295	5,17	108
22	421	3,52	99
24	581	3,16	96

На рис. 1, а–в построены парные диаграммы энергий и сил, действующих на атомы водорода и титана для сравнения MTP-22 и обучающего набора DFT расчетов. Видно, что значения энергий на рис. 1 а группируются в виде обособленных участков. Это объясняется тем, что обучение выполнялось для фиксированного набора концентраций водорода в достаточно узком интервале температур ниже 1000 К. Как следует из рис. 1 б, в, построенные MTP-потенциалы значительно точнее воспроизводят силы, действующие на атомы водорода, в сравнении с силами на атомах титана.

Рис. 1. Сравнение построенного MTP с обучающими данными DFT-расчетов для энергий (а), компонент сил, действующих на атомы водорода (б) и титана (в).

Сравнение результатов MTP с DFT было выполнено для упорядоченного моногидрида γ -TiH и дигидрида TiH 2 . Расчеты проводились без учета влияния конечных температур, выполнялась одноосная деформация вдоль оси c при сохранении объема и условии a = b . На рис. 2 представ-

Физика

лены зависимости потенциальной энергии для различных значений c/a . Решетка моногидрида оказывается стабильной при c>a , с осевым соотношением c/a = 1,10. В то же время для TiH 2 характерен сложный профиль энергии с минимумами в двух точках: c/a = 0,92 и 1,07. Построенный MTP хорошо воспроизводит энергию моногидрида, в то время как для дигидрида энергия оказывается слегка сдвинута относительно DFT-результатов на 0,2 кДж/моль (~2 мэВ/атом), хотя общий вид профиля энергии вдоль деформации c/a предсказывается верно.

Рис. 2. Потенциальная энергия MTP в сравнении с DFT расчетами. для различного тетрагонального искажения c/a кристаллических решеток TiH (а) и TiH 2 (б).

Молекулярно-динамическое моделирование

Расчеты превращения при увеличении температуры для TiH 2 из тетрагональной решетки в кубическую проводились для систем, содержащих 4 000 атомов титана и 8 000 водорода. Для фиксированного значения T система удерживалась в NPT ансамбле на протяжении 40 тыс. шагов молекулярной динамики по 0,5 фс каждый, за вторую половину расчета проходило усреднение длины системы. Внешнее давление баростата задавалось нулевым, причем размеры системы могли изменяться независимо друг от друга, что обеспечило формирование анизотропной тетрагональной структуры. На рис. 3 представлены результаты молекулярно-динамических расчетов, из полученных данных видно, что при температурах выше ~487 K стабильной является кубическая кристаллическая решетка. Полученные результаты несколько выше экспериментальных данных [3, 4], согласно которым ГЦТ-ГЦТ превращение происходит при 310–314 К. Данная разница может объясняться отсутствием в DFT расчётах, использованных для обучения MTP, рассмотрения возбуждения электронных степеней свободы, связанных с конечными температурами.

Рис. 3. МД расчеты с помощью MTP параметров решетки и тетрагонального соотношения c/a дигидрида титана для различных температур

Чирков П.В., Молекулярно-динамическое моделирование стабильности Мирзоев А.А. гидридов титана при помощи машинно-обученных потенциалов

Для гидридов вида TiH x МД расчеты проводились в похожей манере, расчетная система также содержала 4 000 атомов титана при концентрации водорода 1,5 ≤ x ≤ 2,0. Определенное число атомов H, соотвествующее заданному содержанию x , имело случайное распределение по тетраэдрическим междоузлиям ГЦК решетки. Температура в расчетах поддерживалась равной 300 К. Из рис. 4 видно, что образование тетрагональной структуры происходит при содержании водорода x = 1,86, что хорошо согласуется с экспериментальным значением 1,84 [4].

Рис. 4. Параметры решетки и осевое отношение c/a при T =300 K для гидридов с нестехиометрическим составом TiH х .

Заключение

В данной работе исследовались фазовые превращения в системе титан-водород. Для выполнения молекулярно-динамического моделирования на уровне точности первопринципных расчетов было выполнено построение машинно-обученного межатомного потенциала в рамках модели MTP. Для его параметризации в рамках DFT расчетов был создан обучающий набор данных, содержащий 2963 конфигурации. Обученный MTP, содержащий 421 подгоночных коэффициента демонстрирует среднеквадратичное отклонение по энергии 3,52 мэВ/атом и 99 мэВ/Å для проекций сил. Построенный MLIP верифицировался на данных по энергии упорядоченных структур TiH и TiH 2 . Сравнение с первпринципными расчетами показывает хорошее согласие по энергии при одноосной деформации решетки.

Выполненное молекулярно-динамическое моделирование с машинно-обученным межатомным потенциалом в рамках модели MTP позволило оценить температуру ГЦТ-ГЦК превращения для дигидрида TiH 2 , которая составила 487 К, что несколько выше наблюдаемых в экспериментах значений ~310 K [3,4]. Установлено также, что при 300 К тетрагональная решетка образуется для содержаний водорода больших чем 1,86, что согласуется с экспериментальным значением 1,84 [4].

Работа выполнена при поддержке гранта РНФ No. 25-22-20062.