Разработка нового полуэмпирического квантово-химического подхода для моделирования соединений актиноидовНИР

Development of new semiempirical quantum-chemical method for actinide compounds modelling

Источник финансирования НИР

грант РНФ

Этапы НИР

# Сроки Название
1 16 апреля 2019 г.-15 декабря 2019 г. Разработка нового полуэмпирического квантово-химического подхода для моделирования соединений актиноидов
Результаты этапа: Согласно заявленному плану, работу в текущем году проводились по двум основным направлениям: (1) наработке базы данных соединений актиноидов и (2) разработке программного обеспечения для создания полуэмпирической модели. При этом первую задачу можно разбить на две подзадачи: (1.1) непосредственно квантово-химические расчеты, связанные с создание базы данных и (1.2) контроль качества получаемой БД на основе сравнения с существующими в литературе экспериментальными данными Во второй задаче также можно выделить два направления, согласно заявленному набору методов: (2.1) разработка байесовского алгоритма оптимизации полуэмпирического расчетного метода и (2.2) разработка нейросетевого алгоритма оптимизации полуэмпирического расчетного метода. Кроме того, в процессе работы была сформулирована третья задача, ранее предполагаемая на более поздний этап проекта: (3) оптимизация и ускорение существующих (и написанных нами) расчетных алгоритмов. В то же время, относительно большой объем найденных экспериментальных данных позволил на данном этапе использовать исключительно литературные данные для валидации расчетных. 1.1 Наработка базы данных Была создана база данных с рассчитанными на уровне теории связанных кластеров энергиями (dH) более чем 350 реакций с участием актиноидов от актиния до берклия. Однако, для повышения числа и разнообразия возможных структур, работы по расширению базы данных продолжаются. 1.2 Контроль качества БД Поскольку теоретические расчёты проводились для вакуума в качестве среды, то наиболее корректно сравнивать полученные структурные данные с результатами газовой электронографии. Однако, число структур соединений актинидов, исследованных таким образом крайне ограниченно – либо по причинам высокой удельной активности, либо из-за ограниченной летучести соединений. По второй причине единственным более-менее описанным классом соединений являются фториды актинидов. Сопоставленные результаты представлены на рисунке . Для выборки из трёх соединений (ThF4, UF6, NpF6 ) R2 = 0,9989, добавка же UF4 снижает его до 0,9217. Интересным фактом для соединения UF4 является, что рассчитанные значения длин связей U-F не одинаковы: 2,0758; 2,0475; 2,0474; 2,0476 Å. Для трифторидов и тетрафторидов было проведено сравнение с экспериментальными данными о ионных радиусах катионов актинидов. Полученные корреляции представлены на рисунках . Видно, что расчётные данные воспроизводят общие тренды - уменьшение ионного радиуса с увеличением порядкового номера актинида. Из ряда тетрафторидов выбивается уже упомянутый тетрафторид урана – рассчитанное значение оказывается меньше, чем ожидаемое при этом значение ионного радиуса U+4. Для трифторидов из общей тенденции выбивается трифторид америция. 2.1 Байесовский алгоритм Мы адаптировали алгоритм, ранее используемый для настройки гиперпараметров нейронных сетей для оптимизации параметров DFT-функционала. Для этого был написан программный код, состоящий из модулей, ответственных за: чтение исходных структурных и термодинамических данных; расчета энергии (dH) реакции с параметрически-заданным функционалом и экстраполяцией на бесконечный базисный набор; оптимизации параметров функционала и анализа вкладов отдельных элементов функционала. Код опубликован в открытом виде и доступен по адресу github.com/mitrofjr/TPE_func. Код был протестирован на предмет создания функционалов для решения различных общехимических задач. В качестве референсных данных использовалась БД GMTKN55. При этом, тестирование проводилось по двум направлениям: в первом случае рассматривались относительно простые для популярных функционалов задачи, а целью тестирования являлась как общая проверка работоспособности алгоритма, так и оценка вкладов различных функциональных «примитивов» и выявления корреляций между ними для минимизации их числа в итоговом функционале. Во втором случае, рассматривались, наоборот, задачи сложные для популярных DFT функционалов и тестировалась возможность найти оптимальное решение в этом случае. В обоих случаях код справился с задачей. Для вычислительно-простых задач, все отклонения находятся в пределах 1 ккал/моль, что сопоставимо с погрешностью самих референсных данных. Для вычислительно-сложных задач была показана возможность создания функционала, превышающего по точности популярные функционалы (PBE0, B3LYP). Таким образом, была разработана методика и написано программное обеспечения для раработки новых DFT-функционалов под конкретные вычислительные задачи. 2.2. Нейросетевой алгоритм Одним из основных преимуществ нейросетевых алгоритмов, помимо высокой точности предсказания при достаточно большом количестве тренировочных данных, является широкий спектр возможности настройки алгоритма (от выбора архитектуры сети, до подбора параметров обучения) под конкретную задачу. В данном случае было решено разработать наиболее гибкую архитектуру сети, чтобы получить возможность использовать ее как для оптимизации параметров квантово-химических методов, так и для наработки и валидации самой базы данных. Для этой цели мы использовали так называемую «графовую» сверточную нейронную сеть. Рассмотрение графа структуры со свойствами атомов, записанными в узлах графа вместо использования традиционных сетей прямого распространения позволило сделать универсальных инструмент для построения полуэмпирических моделей. Первоначальное тестирование было проведено на относительно больших и популярных наборах данных, где можно было бы сравнить точность алгоритма с аналогами из литературы. Так, при моделировании модуля всестороннего сжатия наш алгоритм показывал близкие значения с лучшим из существующих (R^2 0.91 против 0.93); при моделировании растворимости молекул в воде, наш подход несколько превосходил существующие (ROC AUC 0.97 против 0.93). В случае актинид-содержащих соединений, в качестве исходного набора данных использовалась выборка соединений из Materials project database, включающих в себя актиноиды (порядка 2000 структур). Данная выборка была разбита в соотношении 80/10/10 для дальнейшего использования в качестве обучающей, валидационной и тестовой выборки, соответственно. Итоговое значение средней абсолютной ошибки на тестовой выборки составило 0,017 эВ/атом, что сопоставимо с лучшими предсказательными моделями, обученными на полных базах данных (the Materials Project database, The Open Quantum Materials Database, и т.д.). Таким образом, нами была разработана и протестирована новая и универсальная нейросетевая архитектура, пригодная для моделирования свойств соединений (как молекул, так и кристаллов) актиноидов. 3. Оптимизация существующего ПО Наиболее перпективным вариантом ускорения работы библиотек для квантово-химических расчетов, является их перенос на графические ускорители (GPU) На данный момент удалось внести необходимые изменения в код, используемый для расчета LDA и, частично, GGA-функционалов.. Соответственно, при подключении к квантово-химическим программам новых ускоренных библиотек, расчет данной части взаимодействия осуществляется на графическом процессоре. Таким образом, на данном этапе получилось уменьшить время выполнения части операций, выполняемых при гибридном DFT расчете, что сократит время оптимизации функционала на будущих этапах.
2 1 января 2020 г.-15 декабря 2020 г. Разработка нового полуэмпирического квантово-химического подхода для моделирования соединений актиноидов
Результаты этапа: Был построен новый, оптимизированный для работы с актиноидами, функционал, в основе которого лежит популярный функционал PBE0. Для оптимизации/тестирования использовались реакции, составленные на основе созданной нами базы данных. Абсолютные среднеквадратичные погрешности референсной базы данных: 0.1 А, 12 ккал/моль при воспроизведении геометрии энергии. При этом экспериментальные газофазные данные, взятые для сравнения определены с погрешностями 0.1 А, 10 ккал/моль, соответственно. Таким образом, база данных хорошо воспроизводит существующие экспериментальные данные и подходит в качестве референсной. Релятивистские эффекты учитывались на уровне скалярно-релятивстского гамильтониана Дугласа-Кролла-Гесса и соответствующих полноэлектронных трехэкспоненциальных базисных наборов. Эффекты, связанные с потенциальным многодетерминантным характером волновой функции, на данном этапе не учитывались. Вместо этого, на основе T1-диагностики референсных CCSD(T) расчетов, для создания предварительного функционала были выбраны соединения с пренебрежимо малым многодетерминантным вкладом. Функционал был получен с использованием разработанного нами ранее байесовского подхода при оптимизации параметров включенного градиент-корректировочного функционала PBE, а также долей LDA/GGA корректировки и хартри-фоковского обмена. Параметры полученного функционала приведены в п. 1.3 (иллюстрации в п. 1.6). Погрешности функционала (как по расчету энергий, так и по оптимизации геометрии) согласно кросс-валидации находятся в пределах погрешностей референсной базы данных, в том числе и для заявленных соединений урана. Полученный функционал, согласно первым оценкам, является заметно более точным при работе с соединениями актиноидов, чем существующие. При этом, расчет происходи на порядок быстрее, чем при использовании точных пост-хартри-фоковских-методов. Таким образом, можно заключить, что заявленные научные результаты полностью достигнуты.
3 1 января 2021 г.-15 декабря 2021 г. Разработка нового полуэмпирического квантово-химического подхода для моделирования соединений актиноидов
Результаты этапа: Согласно принятому ранее плану, в 2021 году работы проводились по трем направлениям: 1) Расширение и уточнение существующей базы данных соединений тяжелых элементов, публикация ее в открытом доступе: В ранее созданную в рамках проекта базу данных HERDB (Heavy-Element Reactions DataBase) было добавлено 63 соединения лантаноидов, что существенно расширило её область применения. На данный момент опубликована информация об оптимизированных структурах соединений в газовой фазе, база данных доступна для общего пользования по ссылке https://github.com/SmartChemDesign/HERDB. Аналогично первой версии базы данных, построение включало в себя двухстадийный процесс оптимизации геометрии. Для оптимизации геометрии была использована новая версия программного пакета Orca 5, что позволило решить ряд проблем со сходимостью, в частности для ряда солей органических кислот (рисунок 1). Помимо добавления новых соединений в базу данных, были определены спектральные термы основного состояния для всех новых соединений (таблица 1) и рассчитаны полные внутренние энергии на уровне теории связанных кластеров CCSD(T) для 45 веществ с поправкой на бесконечный базисный набор. Данная часть работы продолжается и на текущий момент. Таблица 1. Полный список соединений, добавленных в базу данных Соединение Терм Соединение Терм LaF3 1S0 La 2D3/2 PrF3 3H4 Ce 3H4 NdF3 4I9/2 Pr 4I9/2 SmF3 6H5/2 Nd 5I4 EuF3 7F0 Pm 6H5/2 GdF3 8S7/2 Sm 7F0 DyF3 6H15/2 Eu 8S7/2 ErF3 4I15/2 Gd 9D2 TmF3 3H6 Tb 6H15/2 LuF3 1S0 Dy 5I8 PrCl3 3H4 Ho 4I15/2 TbCl3 7F6 Er 3H6 LuCl3 1S0 Tm 2F7/2 La(OH)3 1S0 Yb 1S0 Eu(OH)3 7F0 Lu 2D3/2 Gd(OH)3 8S7/2 La3+ 1S0 Tm(OH)3 3H6 Ce3+ 2F5/2 La(NO3)3 1S0 Pr3+ 3H4 Ce(NO3)3 2F5/2 Nd3+ 4I9/2 Pr(NO3)3 3H4 Pm3+ 5H0 Tb(NO3)3 7F6 Sm3+ 6H5/2 Lu(NO3)3 1S0 Eu3+ 7F0 La(HCOO)3 1S0 Gd3+ 8S7/2 Ce(HCOO)3 2F5/2 Tb3+ 7F6 Pr(HCOO)3 3H4 Dy3+ 6H15/2 Nd(HCOO)3 4I9/2 Ho3+ 5I8 Ho(HCOO)3 5I8 Er3+ 4I15/2 Tm(HCOO)3 3H6 Tm3+ 3H6 Lu(HCOO)3 1S0 Yb3+ 2F15/2 Pr(CH3COO)3 3H4 Lu3+ 1S0 Nd(CH3COO)3 4I9/2 Yb(CH3COO)3 2F15/2 Lu(CH3COO)3 1S0 Для первичной оценки необходимости использования многодетерминантного подхода был использован метод NEVPT2, позволяющий учесть как статическую электронную корреляцию за счет учета всех возможных конфигураций активного пространства, так и динамическую корреляцию за счет использования теории малых возмущений Меллера-Плессе. Возможность варьирования размера активного пространства в зависимости от исследуемой системы позволяет в некоторых случаях получать результаты значительно быстрее, чем при использовании метода связанных кластеров. При этом точность учета динамической электронной корреляции достигает сопоставимых величин, как было показано нами ранее. Для нивелирования несовершенства базисного набора была использована схема, включающая оценку ошибки на базе расчетов на уровне теории Меллера-Плессе с базисами Даннинга с различным числом экспонент, описывающих атомные орбитали. Полученная поправка прибавлялась к энергиям, полученным на уровне теории связанных кластеров. Данная схема показывает несколько меньшую точность, чем непосредственное использование теории связанных кластеров, однако последнее практически неприменимо из-за огромных затрат вычислительных ресурсов. Сводные данные о точности полученных подходов представлены в таблице 2. Среднее значение экспериментальной погрешности определения энтальпии реакции составляет, согласно литературным источниками, 7,5 ккал/моль. Расчеты методом связанных кластеров, опубликованные нами ранее, показали погрешность 9 ккал/моль и эти данные были использованы для тренировки DFT функционалов. В рамках данного этапа работ были оценены погрешности при использовании многодетерминантного подхода и экстраполяции на бесконечный базис, составившие 7 и 14 ккал/моль, соответственно. При этом, ухудшение точности при применение схемы экстраполяции может быть объяснено как несовершенством самой схемы, так и высоким качеством базисов SARC, использованных в первой версии базы данных. Таблица 2. Сравнение энтальпий реакций, рассчитанных различными способами с экспериментальными значениями Уравнение реакции CCSD(T)-CBS NEVPT2 CCSD(T) Эксп. UF6 + 2H2O -> UO2F2 + 4HF 47,08 45,92 57,09 45±6 2UO3 + UF6 -> 3UO2F2 -104,95 -60,81 -78,10 -74±12 UF6 + 3H2O -> UO3 + 6HF 123,09 99.29 124,68 104±8 UO3+H2O -> UO2(OH)2 -54,12 -32.82 -48,15 -44±9,5 UO2F2 + 2H2O -> UO2(OH)2 + 2HF 21,89 20.55 19,44 16±2 Таким образом, опубликованная на текущий момент в открытом доступе версия базы данных включает информацию о 176 атомах, ионах и солях лантаноидов и актиноидов. Для каждого соединения представлены данные о геометрии, рассчитанной на уровне теории Меллера-Плессе (MP2) с учетом скалярных релятивистских поправок, заряде и электронном терме. Для большинства соединений так же присутствуют данные о полной внутренней энергии и электронной структуре (в виде molden-файла, содержащего рассчитанные орбитали) рассчитанные на уровне теории связанных кластеров CCSD(T). База постепенно пополняется также описанными выше данными об электронной структуре и полной внутренней энергии, с учетом несовершенства базисного набора. В настоящее время готовится к публикации статья в рецензируемом журнале, посвященная второй версия базы данных. 2) Систематизация и публикация результатов по разработанному методу генерации конформаций и конформационного анализа: Ранее, для автоматизации процесса поиска глобального минимума на поверхности потенциальной энергии конформаций комплексов f-элементов нами была предложена схема, основанная на стохастической генерации конформаций и их дальнейшей двухэтапной оптимизации, состоящей из оптимизации конформаций самого лиганда и его координации относительно иона металла. Для относительно малых и негибких неорганических анионов, ожидаемо более важным был второй этап, позволивший существенно автоматизировать и ускорить процесс получения структур соединений для базы данных. Для более гибких органических лигандов, в свою очередь, гораздо важнее оказался первый этап предорганизации самого лиганда. Совмещений стохастической генерации конформаций на базе матрицы связности и их локальной предоптимизации позволило не только находить их равновесную геометрии, но и получать общий вид их поверхности потенциальной энергии (ППЭ), что является даже более важной задачей ввиду возможных различий в связывающей конформации лиганда и его конформацией в глобальном минимуме ППЭ. Так, например, молекула бипиридина, связывающая металл двумя азотами, находящимися в син-конформации, в свободном состоянии склонна принимать транс-конформацию, непригодную для образования комплекса. Сама важность конформационной подвижности лиганда для его комплексообразования и необходимость оценки высоты энергетических барьеров ППЭ обсуждалась в литературе в рамках моделирования комплексообразования и ранее. Основным подходом к моделированию являлась последовательная оптимизация геометрии соединения при изменении с постоянным шагом величины двугранных углов, соответствующих подвижным (например, одинарным C-C) связям в молекуле. При этом, среди основных сложностей подобного моделирования можно отметить высокое вычислительное время, увеличивающееся по степенной функции от числа степеней свободы; зависимость от выбранной пользователем начальной геометрии и возможность системы прийти в состояние «храповика», когда последовательное вращение и оптимизация геометрии приводят систему в состояние, выход из которого не может быть осуществлен путем продолжения вращения и оптимизации. Разработанный нами ранее подход лишен указанных недостатков, а рост вычислительных требований линеен относительно числа степеней свободы, что делает подобные расчеты существенно короче. При использовании метода для пополнения базы данных комплексами f-элементов и малых органических молекул, было отмечено, что, для более гибких лигандов, метод находит и связывающие, и свободные конформации. Причем, для разных лигандов, были найдены отличающиеся значения разности внутренних энергий для указанных конформаций, коррелирующие с известными константами устойчивости комплексов. Сам расчет внутренней энергии проводился на относительно простом уровне теории функционала плотности с гибридным функционалом PBE0 и поправкой D3 для учета дисперсионных взаимодействий. Таким образом, метод показал перспективность оценки эффективности комплексообразования без необходимости использования вычислительно дорогих расчетов комплексов. Для подтверждения гипотезы была рассмотрена серия лигандов, на основе фосфин-оксидных и моноамидных производных пиридина, отличающихся только заместителями (рис. 2). С использованием разработанного метода конформационного анализа, были получены структуры и соответствующие значения полной внутренней энергии лигандов в связывающем и свободном состоянии. Методом 2D-ЯМР (NOESY) было показано, что метод корректно находит глобальный минимум на поверхности потенциальной энергии лиганда. Для учета возможного влияния растворителя также были рассчитаны энергии сольватации каждой из конформаций и предложена сумма энергий сольватации и предорганизации как новый параметр оценки эффективности комплексообразования (SCPE – solvent-corrected preorganization energy). Экспериментальная проверка показала хорошую корреляцию полученного параметра с значениями эффективности связывания Eu(III) и Am(III) (рис. 3) Коэффициенты корреляции R^2 составили 0.8-0.85 для случая включения или не включения в рассмотрение несимметричных лигандов. При этом, полученные данные были использованы для анализа связи структуры лиганда и эффективности связывания металла, а полученный алгоритм имплементирован в процесс создания базы данных. Результаты были опубликованы в статье в Journal of molecular liquids [10.1016/j.molliq.2020.115098] 3) Доработка полученного функционала: Для получения итоговых функционалов были использованы три базовых существующих функционала: PBE0, рассмотренный нами ранее, B3LYP, являющийся одним из наиболее популярных и часто используемых функционалов в том числе для соединений f-элементов и TPSSh, показавший хорошие результаты на предыдущих этапах исследования. Для улучшения работоспособности функционалов также использовались две схемы – добавление учета электронной корреляции методом малых возмущений в корреляционную часть обменно-кореляционного функционала (создание т.н. двойного гибридного функционала) и разделение Хартри-Фоковского обмена на отдельные участки, отвечающие за локальные и нелокальные взаимодействия (создание т.н. range-separated функционала). При этом, одновременно оптимизировались как внутренние параметры функционала, так и параметры его модернизации, позволяя оценить оба пути в рамках одного тренировочного запуска. Оптимизация проводилась с использованием опубликованного нами ранее алгоритма [10.1021/acs.jpca.9b09093] модифицированного для использования программного пакета Orca для расчета полных внутренних энергий на каждом этапе. Тренировка проводилась в течение 600-800 итераций для различных подходов, в зависимости от требуемого на расчет времени. В наилучшем случае, для достижения уровня погрешности, соответствующего погрешности исходной базы данных, потребовалось порядка 500 итераций. Каждой итерации соответствовал расчет энтальпии трех построенных на основе базы данных реакций. За счет стохастической выборки тренировочных наборов достигалась частичная внутренняя валидация результатов для предотвращения переобучения на тренировочных данных. Итоговое тестирование проводилось на существующих экспериментальных данных. Ход тренировки (т.н. кривые обучение), анализ вкладов параметров и их корреляции показан на рис. 4. Из протестированных подходов наихудшие результаты показало использование range-separated схемы. Наилучшая комбинация, полученная в процессе тренировки данной модификации функционала, с высокой точностью соответствовала уже существующему CAM-B3LYP функционалу, что позволило нам исключить этот путь из дальнейших рассуждений. С другой стороны, из исследованных вариантов тренировки, модифицированной версии функционала PBE0 удалось пройти рубеж погрешности референсной базы данных в отношении моделирования изменения энтальпии реакций. При этом, дополнительное тестирование функционала на выборке реакций W4-11 из базы данных GMTKN55 не выявило существенных отличий в точности с оригинальным функционалом PBE0, что позволяет предполагать широкую область применения функционала. Результаты сравнения точности полученных функционалов с существующими относительно моделирования реакций актиноидов представлены на рис. 6. По итогам тренировки, был получен вариант гибридного функционала, позволяющий проводить расчеты энергии с точностью, сравнимой с референсной точностью базы данных (порядка 9 ккал/моль). Более того, использование более вычислительно-тяжелых схем не позволило существенно улучшить точность расчетов. Дальнейшее тестирование данного функционала проводилось в режиме оптимизации геометрий существующих в газовой фазе соединений актиноидов. На рис. 7 видно, что разработанный подход показал в целом удовлетворительные значения воспроизведения геометрий соединений в целом и длин связей актинид-неметалл в частности, средние среди исследованных методов. Помимо разработки нового функционала, в процессе исследования была показана возможность использования существующего функционала TPSSh, точно воспроизводящего экспериментальные геометрии, и уступающему в оценки энтальпии реакций. Таким образом, в рамках исследования были рассмотрены различные подходы к формированию обменно-корреляционного функционала. В итоге был получен новый гибридный функционал relPBE0, который не привносит дополнительную ошибку в расчет полной внутренней энергии по отношению к базовым CCSD(T) расчетам. Результаты опубликованы в журнале Journal of Chemical Physics [10.1063/5.0067631]
4 1 января 2022 г.-15 декабря 2022 г. Разработка нового полуэмпирического квантово-химического подхода для моделирования соединений актиноидов
Результаты этапа: 1) Вторая версия базы данных, включающая оценку внутренней энергии различными методами, уточненное тестирование полученных значений относительно известных газофазных реакций и структур как соединений актиноидов, так и лантаноидов. 2) Вторая версия функционала, включающая его тренировку и тестирование на данных как о реакциях актиноидов, так и лантаноидов а также расширенную информацию о привносимых им погрешностях, включая поведение с различными базисными наборами

Прикрепленные к НИР результаты

Для прикрепления результата сначала выберете тип результата (статьи, книги, ...). После чего введите несколько символов в поле поиска прикрепляемого результата, затем выберете один из предложенных и нажмите кнопку "Добавить".