Tetraoxa[8]circulene monolayer as hydrogen storage material: model with Boys–Bernardi corrections within density functional theory
- Authors: Anikina E.V.1, Babailova D.V.1, Zhilin M.S.1, Beskachko V.P.1
-
Affiliations:
- Institute of Natural Sciences and Mathematics, South Ural State University
- Issue: No 1 (2024)
- Pages: 23-32
- Section: Articles
- URL: https://ogarev-online.ru/1028-0960/article/view/256986
- DOI: https://doi.org/10.31857/S1028096024010049
- EDN: https://elibrary.ru/DPYRHN
- ID: 256986
Cite item
Full Text
Abstract
The parameters of molecular hydrogen adsorption on a tetraoxa[8]circulene monolayer were studied using the density functional theory with dispersion interaction corrections (semi-empirical and analytical). The calculations were carried out using two different approaches to the system wave function representation: atomic-like orbital basis set and plane wave basis. Utilizing a less computationally expensive pseudoatomic basis, it is possible to obtain results for molecular hydrogen adsorption consistent with values calculated with plane waves if the atomic-like basis is optimized and basis set superposition error is corrected for both hydrogen binding energy and geometrical characteristics. Otherwise, the H2 binding energy will be overestimated by 4–6 times (sometimes even more, by 20); and the hydrogen-monolayer distance will be underestimated by 10-20%. The obtained optimized parameters of the pseudoatomic basis set can be used for further study of the modified forms of the tetraoxa[8]circulene monolayer. Moreover, our calculations showed that the hydrogen binding to a pristine tetraoxa[8]circulene monolayer is predominantly van der Waals with an energy of 60–90 meV, which is several times less than the desired range of 200–600 meV. To achieve such values, it will be necessary to modify the surface of the monolayer, creating more active sorption cites, for example, by decorating it with metals or applying structural defects.
Full Text
ВВЕДЕНИЕ
Материалы, размерность которых можно считать пониженной (меньше 3), нередко проявляют уникальные свойства благодаря своей предельно высокой удельной поверхности. Если такие материалы имеют еще и низкую плотность, а также химическую и тепловую стойкость, то они являются предметом пристального внимания как материалы для хранилищ водорода [1]. Типичным примером являются низкоразмерные аллотропы углерода [2]. В современных базах кристаллографических данных присутствуют сотни подобных материалов, и этот список непрерывно пополняется новыми кандидатами. К тому же, поверхность таких материалов можно модифицировать разными способами. Соответственно, экспериментальная проверка пригодности таких материалов для хранения водорода достаточно трудоемкая и дорогостоящая. Поэтому на этапе поиска структуры с заданными свойствами эффективным является применение компьютерного моделирования, которое позволяет относительно быстро протестировать множество вариантов и отобрать для экспериментального получения только те материалы, которые обладают наиболее подходящими характеристиками. Без предварительных эмпирических данных исследовать новые системы позволит теория функционала электронной плотности, но применение конкретных программных пакетов, в которых эта теория имплементирована, сопряжено с выбором приближений и параметров моделирования, которые могут сильно влиять на точность и быстроту получения результатов.
Ключевым приближением, напрямую влияющим как на скорость расчета, так и на получаемый результат, является выбор вида конечного базиса, в котором будет представлена волновая функция исследуемой системы. Быстрые результаты даже при условии ограниченных вычислительных ресурсов (например, только при наличии ПК) можно получить с помощью базиса локализованных атомноподобных орбиталей (АО) [3]. Для этого базиса, в отличие от базиса плоских волн (ПВ), не является проблемой наличие в модели больших объемов вакуума, неизбежно возникающих при исследовании сорбции с использованием периодических граничных условий, поскольку базисные орбитали хотя и протяжены, но ограничены в пространстве. С другой стороны, эта ограниченность в совокупности с зависимостью АО-базиса от моделируемой атомной структуры объекта приводят к источнику погрешности метода, известному как ошибка суперпозиции базисного набора (BSSE).
Эта ошибка возникает всякий раз, когда требуется вычислить энергию связи между отдельными частями моделируемой системы (например, энергию связи молекулы H2 с материалом основы). Когда исследуемые подсистемы находятся далеко друг от друга (не взаимодействуют), то каждая из них будет описана только “своими” орбиталями, но как только подсистемы сближаются, их базисные функции начинают перекрываться, и тогда в описании одной подсистемы будут участвовать орбитали второй подсистемы. Такое “излишнее” описание приводит к искусственному притяжению между подсистемами, что приведет и к завышенным величинам энергии связи, и к слишком близкому положению подсистем после стандартной оптимизации геометрии системы до достижения порога силы. В слабосвязанных структурах на фоне небольшой энергии связи (не более 100 мэВ/H2) ошибка суперпозиции базисного набора может быть значительной [4], и заметное влияние она может оказывать и на исследования объемных кристаллов [5].
Монослой тетраоксо[8]циркулена (ТОС), выбранный нами в качестве сорбента молекулярного водорода, имеет пористую структуру и термически устойчив, как показывают расчеты [6]. Он относится к новому классу органических полупроводников и впервые был синтезирован лишь недавно [7]. Первоначальный интерес к нему возник в связи с возможными применениями в области гибкой электроники и квантовой информатики (как возможному хранилищу кубитов) [6, 8, 9], а попыток оценить его как материал для хранилищ водорода нам обнаружить не удалось. Таким образом, цель настоящей работы заключалась в оценке свойств нового материала как материала водородных хранилищ и в разработке методики повышения надежности этой оценки с использованием только быстрых методов моделирования.
МОДЕЛИ И ПАРАМЕТРЫ МОДЕЛИРОВАНИЯ
Для расчетов в рамках теории функционала электронной плотности использовали программные пакеты SIESTA [10] и VASP [11]. В первом волновая функция системы представлена в базисе локализованных атомноподобных орбиталей, генерируемых по схеме Санки (Sankey) [12], во втором – в базисе плоских волн. Спин-поляризованные расчеты как чистого, так с адсорбированным водородом монослоя TOC в пакете VASP показали, что магнитный момент системы равен 0, поэтому в пакете SIESTA для экономии компьютерных ресурсов проводили расчеты без спиновой поляризации. Для учета дисперсионного взаимодействия водорода с монослоем, которое может играть заметную роль в связывании молекул водорода с низкоразмерным материалом на основе углерода [13], использовали как полуэмпирические поправки Гримме второго (D2 [14]) и третьего (D3 [15]) поколения к обменно-корреляционному функционалу (ОКФ) Пердью–Берка–Эрнзерхофа (Perdew–Burke–Ernzerhof, PBE) [16], так и аналитический ван-дер-ваальсовый ОКФ Берланда–Хильдгаарда (Berland-Hyldgaard, vdW-BH) [17, 18]. Все расчеты проводили в приближении псевдопотенциала, для моделирования в пакете SIESTA псевдопотенциалы взяты из базы Fritz-Haber-Institute [19], в пакете VASP использовали фирменные псевдопотенциалы версии 2012 года. Энергия обрезки базиса плоских волн в пакете VASP составляла 600 эВ. В качестве валентных рассматривали 2s22p2 электроны для атомов углерода, 2s22p4 – для атомов кислорода и 1s – для водорода. Оптимизированные параметры моделирования [20], характеризующие схему вычислений, представлены в табл. 1. Они гарантировано позволяют получить погрешность вычислений энергии связи молекулы водорода порядка 1 мэВ.
Таблица 1. Параметры моделирования
Программный пакет | SIESTA | VASP |
Приближение для ОКФ | PBE+D2; vdW-BH | PBE+D2; PBE+D3 |
Разбиение обратного пространства, k-точки | 9×9×1 | 7×7×1 |
Разбиение пространства: параметры MeshCutoff (SIESTА) и PREC (VASP) | 400 Рб (5.44 кэВ) | Accurate |
Критерий сходимости полной энергии, эВ | 10-6 | |
Критерий сходимости силы, эВ/Å | 1.29×10–3 | 10–3 |
Размеры ячейки (в скобках – значения для базиса по умолчанию), Å3 | 8.41×8.41×40 (8.40×8.40×40) | 8.39×8.39×25 |
Тетрагональная ячейка тетраоксо[8]циркулена содержала 24 атома: 20 атомов углерода и 4 атома кислорода (рис. 1). Из-за использования периодических граничных условий потребовалась достаточно большая вакуумная прослойка в направлении, перпендикулярном монослою, чтобы можно было пренебречь взаимодействием структуры с ее изображением. Оптимальный параметр трансляции в плоскости монослоя вычисляли по минимуму полной энергии системы. Размеры ячейки для разных ОКФ и программных пакетов приведены в табл. 1. Полученные в настоящей работе параметры трансляции согласуются с предыдущими исследованиями TOC [6].
Рис. 1. Атомная структура ячейки моделирования монослоя тетраоксо[8]циркулена. Атомы углерода обозначены серым цветом, атомы кислорода – красным. Границы ячейки моделирования отмечены черной сплошной линией. Все изображения структур получены с помощью программного пакета VESTA 3. Длины связей после оптимизации геометрии в пакете VASP (в скобках – значения после оптимизации в SIESTA): 1 – 1.38 (1.39); 2 – 1.41 (1.42); 3 – 1.40 (1.40); 4 – 1.42 (1.43); 5 – 1.41 (1.41); 6 – 1.46 (1.46) Å.
Энергию связи водорода с тетраоксо[8]циркуленом рассчитывали следующим образом:
Ebind = ETOC + EH2 – ETOC + H2 – |ECP|, (1)
где ETOC + H2 – полная энергия монослоя с адсорбированной молекулой водорода; ETOC – полная энергия изолированного монослоя тетраоксо[8]циркулена; – полная энергия изолированной молекулы водорода; ECP – поправка Бойса–Бернарди [21] к ошибке суперпозиции базисного набора (в ПВ-базисе равна 0). Заметим, что положительное значение энергии связи отвечает притяжению молекулы водорода к монослою. Поправку ECP рассчитывали как в [22].
РЕЗУЛЬТАТЫ
Оптимизация атомноподобного базиса
В пакете SIESTA АО-базис вида DZP (double-ζ polarized) [23] оптимизировали по процедуре, аналогичной описанной в работе [22], с тем отличием, что в качестве геометрического критерия оптимальности использовали “давление” P, зависящее от разницы между текущим параметром трансляции и его равновесным значением. Отрицательные значения P свидетельствуют о том, что размер ячейки моделирования необходимо уменьшить, положительные – увеличить. Для водорода использовали следующие параметры: радиус обрезки r = 6.05a0 (a0 – радиус Бора, 0.529 Å), радиус модифицированной орбитали rm = 2.52a0. Зависимости характеристик TOC от параметров орбитали O2p представлены на рис. 2. Как видно из рис. 2, при увеличении основного радиуса обрезки r параметры системы (здесь полная энергия и заряд, перешедший на атомы кислорода, рассчитанный по процедуре Милликена) устанавливаются около определенного уровня. Поэтому в качестве оптимального выбирали минимальное значение r, при котором колебания параметров незначительны, что обеспечит некоторый баланс между точностью вычислений и затрачиваемыми вычислительными ресурсами. Радиус модифицированной орбитали rm, задаваемый с помощью параметра PAO.SplitNorm [23], выбирали по минимуму полной энергии системы. Эту процедуру применяли ко всем орбиталям, рассматриваемым в данной АО-модели тетраоксо[8]циркулена. Полученные таким способом оптимизированные параметры АО представлены в табл. 2 вместе с параметрами, устанавливаемыми программой SIESTA по умолчанию. Из табл. 2 видно, что радиусы обрезки r, устанавливаемые по умолчанию, во всех случаях гораздо меньше оптимизированных радиусов (для кислорода, например, более чем в два раза), что приводит к значительной ошибке суперпозиции базисного набора, в чем мы убедимся ниже.
Рис. 2. Зависимости полной энергии ячейки монослоя тетраоксо[8]циркулена Etot, параметра P и заряда, перешедшего на атом кислорода QO от: (а) радиуса обрезки r орбитали О2p; (б) параметра SplitNorm (задает радиус обрезки rm модифицированной орбитали). Заряд кислорода приведен у основания столбцов (1). a0 = 0.529 Å.
Таблица 2. Параметры базисного набора для орбиталей C2s, C2p, O2s и O2p. a0 = 0.529 Å. В скобках приведены значения, генерируемые для базиса DZP пакетом SIESTA по умолчанию
Орбиталь | r, a0 | SplitNorm | rm, a0 |
C2s | 7.09 (4.09) | 0.35 (0.15) | 2.81 (3.35) |
C2p | 6.57 (4.87) | 0.25 (0.15) | 3.18 (3.48) |
O2s | 7.00 (3.31) | 0.15 (0.15) | 2.57 (2.51) |
O2p | 8.13 (3.94) | 0.20 (0.15) | 2.48 (2.54) |
Влияние ошибки суперпозиции базисного набора на геометрию молекулы водорода
Перейдем к моделированию сорбции водорода. Для начала посмотрим, как характеристики атомноподобного базисного набора влияют на сорбцию H2. Ранее уже было показано [22], что энергия связи водорода с основой может быть значительно переоценена из-за ошибки суперпозиции базисного набора, и для получения корректного значения в случае слабосвязанных систем необходимо учитывать поправку к этой ошибке. Однако из-за переоцененного притяжения между подсистемами моделируемой структуры (в данном случае – молекулы водорода и монослоя TOC) также получатся и некорректные геометрические характеристики системы. Например, при стандартной процедуре оптимизации геометрии (до достижения порога точности по силе) можно ожидать сокращения расстояния между молекулой H2 и монослоем.
Рассмотрим подробно моделирование сорбции водорода в АО-базисе с параметрами, установленными по умолчанию, на примере одной конфигурации системы, когда молекула H2 располагается над порой, образованной атомами кислорода (конфигурация О1) в приближении Пердью–Берка–Эрнзерхофа с учетом поправки Гримме второго поколения (PBE+D2). Сначала проводили стандартную оптимизацию геометрии по силовому критерию (табл. 1). Затем полученную конфигурацию всех атомов TOC фиксировали и вручную изменяли только расстояние водород–TOC с шагом 0.05 Å как в большую, так и в меньшую сторону от полученного после оптимизации значения. На каждом шаге не проводили какой-либо релаксации структуры (иначе бы атомы вернулись к исходной конфигурации) и вычисляли энергию связи водорода по формуле (1) для двух случаев: с учетом и без учета поправки Бойса–Бернарди (для построения зависимости без учета поправки считали ECP = 0 эВ). Полученные зависимости энергии связи от расстояния d между центром молекулы H2 и монослоем TOC без и с учетом поправки Бойса–Бернарди представлены на рис. 3. Расстояние, соответствующее минимуму зависимости с рис. 3, и будет равновесным расстоянием для использованного метода подсчета энергии связи (с учетом поправки и без). Видно, что без учета поправки к ошибке суперпозиции базисного набора расстояние d получилось на 0.6 Å меньше, что составляет более 30% от этого значения d. Очевидно, что нельзя игнорировать ошибку суперпозиции базисного набора и при расчете геометрических характеристик слабосвязанных систем. По зависимости Ebind(d) можно также оценить погрешность, с которой определено расстояние d после учета поправки Бойса–Бернарди. Если погрешность вычисления энергии связи ΔEbind ~ 1 мэВ, тогда такое изменение энергии произойдет в Ebind(d) на расстоянии Δd от минимума, равном примерно 0.1 Å для данных, показанных на рис. 3.
Рис. 3. Зависимости энергии связи водорода в конфигурации О1 от расстояния d между центром молекулы H2 и монослоем TOC, вычисленные до и после внесения поправки к ошибке суперпозиции базисного набора в геометрию структуры. Результаты приведены для приближения PBE+D2 и атомноподобного базиса, генерируемого по умолчанию в пакете SIESTA.
Осталось понять, какой из способов получения равновесного расстояния d дает корректный результат и как на это влияет процедура оптимизация базиса и выбор обменно-корреляционного функционала. Для выяснения этого вопроса далее представлены результаты моделирования сорбции водорода в пакете SIESTA и их сравнение с результатами расчетов в ПВ-базисе, полученными в пакете VASP, для которых поправка Бойса–Бернарди равна 0.
Влияние оптимизации атомноподобного базисного набора на ошибку суперпозиции базисного набора
Сначала наиболее устойчивые конфигурации молекул водорода были получены в ПВ-базисе. Соответствующие структуры, смоделированные в приближении PBE+D2, представлены на рис. 4. Затем эти конфигурации были использованы в качестве стартовых для расчетов в пакете SIESTA (с приближениями PBE+D2 и vdW-BH для ОКФ). После оптимизации геометрии в приближении PBE+D2 сохранились все рассмотренные конфигурации молекул водорода при использовании и базиса по умолчанию, и оптимизированного базиса (менялась только величина d). В приближении vdW-BH конфигурация С1⊥ перешла в С1||, а в базисе по умолчанию, кроме этого, конфигурация О2 перешла в О1. Вычисленные энергии связи и расстояния между монослоем TOC и водородом в пакете SIESTA приведены в табл. 3. Для всех приведенных конфигураций погрешность вычисления расстояния после учета поправки Бойса–Бернарди к геометрии составила Δd ≈ 0.1 Å.
Рис. 4. Оптимизированные в приближении PBE+D2 в ПВ-базисе атомные структуры тетраоксо[8]циркулена с адсорбированной молекулой водорода. Атомы углерода, кислорода и водорода обозначены серым, красным и зеленым цветом соответственно. Все расстояния d между монослоем и центром молекулы H2 приведены в Å.
Таблица 3. Энергетические и геометрические характеристики молекулы H2, адсорбированной на TOC, рассчитанные в пакете SIESTA. Здесь ΔdCP – увеличение расстояние между водородом и TOC после поправки в геометрии
ОКФ | Конфигурация | До поправки в геометрии | После поправки | ||||
Ebind, мэВ | d, Å | |ECP|, мэВ | Ebind, мэВ | d, Å | |ECP|, мэВ | ||
PBE+D2 | Базис по умолчанию | ||||||
С1|| | 60 | 2.36 | 231 | 70 | 0.25 | 201 | |
С1⊥ | 65 | 2.58 | 178 | 69 | 0.17 | 166 | |
O1 | 28 | 1.76 | 561 | 101 | 0.60 | 352 | |
C2 | 54 | 2.54 | 218 | 66 | 0.25 | 184 | |
O2 | 55 | 2.50 | 205 | 67 | 0.23 | 171 | |
Оптимизированный базис | |||||||
С1|| | 68 | 2.54 | 46 | 71 | 0.15 | 39 | |
С1⊥ | 51 | 2.73 | 36 | 53 | 0.13 | 32 | |
O1 | 74 | 2.26 | 16 | 75 | 0.10 | 14 | |
C2 | 63 | 2.72 | 47 | 66 | 0.15 | 39 | |
O2 | 65 | 2.68 | 45 | 68 | 0.12 | 39 | |
vdW-BH | Базис по умолчанию | ||||||
С1|| | 45 | 2.53 | 222 | 80 | 0.65 | 106 | |
O1 | 26 | 1.90 | 539 | 136 | 0.75 | 258 | |
C2 | 38 | 2.68 | 209 | 71 | 0.60 | 104 | |
Оптимизированный базис | |||||||
С1|| | 58 | 2.77 | 36 | 65 | 0.28 | 20 | |
O1 | 83 | 2.73 | 8 | 84 | 0.07 | 8 | |
C2 | 54 | 2.91 | 38 | 61 | 0.27 | 22 | |
O2 | 59 | 2.90 | 36 | 64 | 0.23 | 25 |
Из табл. 3 видно, что, во-первых, влияние ошибки суперпозиции базисного набора на величину энергии связи очень большое: до оптимизации базиса поправка Бойса–Бернарди составляет 75–85% и может доходить до 95% от величины нескорректированной Ebind. Ситуация немного улучшается после учета ошибки суперпозиции базисного набора в геометрии структуры (особенно в приближении vdW-BH, где доля ошибки падает до ~60%), однако значительно (в ~2 раза) уменьшить долю поправки можно только при использовании оптимизированного базиса. Наиболее чувствительной к оптимизации базиса оказалась конфигурация О1, где в базисе по умолчанию величина |ECP| превышала 0.5 эВ, а в оптимизированном базисе, наоборот, она была минимальной из полученных. Во-вторых, после оптимизации базиса гораздо меньшее влияние ошибка суперпозиции базисного набора оказывает на d: даже до учета поправки в геометрии во всех конфигурациях d получилось заметно больше, чем в базисе по умолчанию, а относительное изменение ΔdCP/d после поправки в геометрию падает в ~2 раза по отношению к таковому в базисе по умолчанию. Вновь наиболее чувствительной к влиянию ошибки суперпозиции базисного набора на геометрию и оптимизации базиса оказалась конфигурация О1: в базисе по умолчанию у нее максимальное значение ΔdCP, а в оптимальном – минимальное из полученных. Интересно, что именно эта конфигурация оказалась наиболее выгодной энергетически в обоих приближениях ОКФ после поправки к ошибке суперпозиции базисного набора в геометрии. Следует отметить, что до оптимизации базиса и до внесения поправки в d оба использованных приближения для ОКФ предсказывали неверную наиболее выгодную конфигурацию водорода на монослое TOC. Характеристики, вычисленные с помощью разных функционалов, качественно похожи, только в приближении vdW-BH получается заметно большее расстояние между H2 и TOC, и в этом же приближении чувствительность параметров сорбции водорода к поправке ошибки суперпозиции базисного набора в геометрию структуры выражена сильнее всего.
Таким образом, очевидно, что в слабосвязанных системах ошибка суперпозиции базисного набора играет ключевую роль при оценке как энергетических, так и геометрических характеристик структуры. Но дает ли поправка Бойса–Бернарди более точные результаты – отдельный вопрос. Чтобы прояснить это, сравним результаты моделирования в пакете SIESTA с таковыми в пакете VASP, где, как уже упомянуто выше, поправка Бойса–Бернарди равна 0. В табл. 4 представлены результаты расчетов, сделанных в приближении PBE+D2 (этот ОКФ реализован в обоих пакетах).
Таблица 4. Энергетические и геометрические характеристики адсорбированной молекулы водорода в приближении PBE+D2 в зависимости от использованного базисного набора. Здесь ∆Ebind = – , ε = (dSIESTA – dVASP)/dVASP
Базис | Атомноподобные орбитали по умолчанию | Оптимизированные атомноподобные орбитали | Плоские волны | |||||||
До поправки в геометрии | После поправки | До поправки в геометрии | После поправки | |||||||
Конфигурация | ΔEbind, мэВ | ε, % | ΔEbind, мэВ | ε, % | ΔEbind, мэВ | ε, % | ΔEbind, мэВ | ε, % | ΔEbind, мэВ | d, Å |
С1|| | -2 | -11.9 | 8 | -2.6 | 6 | -5.2 | 9 | 0.4 | 62 | 2.68 |
С1⊥ | 8 | -8.5 | 12 | -2.5 | -6 | -3.2 | -4 | 1.4 | 57 | 2.82 |
O1 | -46 | -25.1 | 27 | 0.4 | 0 | -3.8 | 1 | 0.4 | 74 | 2.35 |
C2 | -9 | -12.4 | 3 | -3.8 | 0 | -6.2 | 3 | -1.0 | 63 | 2.90 |
O2 | -8 | -10.7 | 4 | -2.5 | 2 | -4.3 | 5 | 0.0 | 63 | 2.80 |
Видно, что до оптимизации АО-базиса есть риск получить результаты, слишком далекие от значений в ПВ-базисе: для конфигурации О1 отклонения как в энергии связи, так и в расстоянии H2–TOC (до поправки в геометрию) слишком велики, составляют до ~60% от ФОРМУЛА и более 25% от dVASP. Однако даже в базисе по умолчанию после учета поправки в геометрию получились расстояния, близкие к результатам VASP (относительная разница ε упала в более, чем 3 раза, и не превышала 4%). После оптимизации базиса адекватные энергии связи можно получить и до внесения поправки в геометрию. Полного совпадения результатов ожидать не стоит, так как в пакетах SIESTA и VASP использованы разные псевдопотенциалы, не говоря уже о схемах и алгоритмах вычислений и других деталях. После учета поправки в d отличия с геометрией, полученной в ПВ-базисе составили ~1%. Таким образом, при использовании АО-базиса следует: с осторожностью использовать короткие орбитали и лучше предварительно оптимизировать параметры этого базиса; обязательно учитывать поправку Бойса–Бернарди при вычислении энергии связи; для получения геометрии, хорошо согласующейся с ПВ-расчетами, необходимо и вычисление расстояния d проводить с учетом поправки к ошибке суперпозиции базисного набора.
Сорбция водорода на монослое тетраоксо[8]циркулена
Сравним теперь результаты, полученные с помощью различных поправок на дисперсионное взаимодействие. Поскольку в приближении PBE+D2 величины энергий связи и расстояний, вычисленные с помощью АО- и ПВ-базисов, практически совпадают после оптимизации базиса и учета поправки Бойса–Бернарди в энергии и геометрии, для дальнейшего анализа мы использовали результаты расчетов в пакете VASP. Итоговые результаты представлены на рис. 5. Конфигурация С1⊥ была исключена из рассмотрения, т.к. в расчетах с полуэмпирическими поправками она имела минимальную энергию связи, а при использовании аналитического функционала оказалась неустойчивой и переходила в С1||. На рис. 5 обозначение С1 соответствует конфигурации С1||.
Рис. 5. Энергия связи Ebind и расстояние между центром H2 и монослоем TOC, вычисленные с дисперсионными поправками для четырех наиболее энергетически выгодных конфигураций. Результаты с поправками Гримме получены в пакете VASP, с аналитическим vdW функционалом – в пакете SIESTA.
Из рис. 5 видно, что все дисперсионные поправки дают качественно похожие результаты. Максимальные оценки энергии связи получились при использовании полуэмпирических поправок Гримме третьего поколения, для некоторых конфигураций (С1 и О1) эти оценки оказались примерно на 20% больше, чем при использовании других приближений. Такие результаты согласуются с предыдущими исследованиями углеродосодержащих низкоразмерных наноматериалов [24]. Аналитический ван-дер-ваальсов функционал дал близкие к PBE+D2 энергии связи и самые большие расстояния d. Наиболее компактное расположение водорода наблюдали при использовании поправок Гримме второго поколения, хотя энергии связи с этими поправками, как правило, были минимальными. Как и ожидалось, энергетически наиболее выгодное расположение водорода оказалась в поре, образованной атомами кислорода (О1). Но даже в таком положении энергия связи не превышала ~90 мэВ, что слишком мало, чтобы обеспечить эффективные циклы сорбции/десорбции водорода в хранилищах при температурах близких к комнатной. Поэтому на следующем этапе исследований необходимо попытаться создать на поверхности монослоя тетраоксо[8]циркулена более активные центры сорбции, например, легируя TOC атомами металлов. Для ряда легких металлов уже было показано, что они хорошо связываются с монослоем [6].
ЗАКЛЮЧЕНИЕ
В настоящей работе с помощью методов теории функционала электронной плотности исследована сорбция молекулярного водорода на монослое тетраоксо[8]циркулена. Расчеты сделаны с использованием базисов локализованных атомно-подобных орбиталей (пакет SIESTA) и плоских волн (пакет VASP) c учетом поправок на дисперсионное взаимодействие разного вида (полуэмпирических и аналитических).
Показано, что для получения адекватных энергетических и геометрических характеристик слабосвязанных систем в АО-базисе и согласования этих результатов с полученными в ПВ-базисе необходима оптимизация параметров АО-базиса, а также учет поправки Бойса–Бернарди к ошибке суперпозиции базисного набора при вычислении как энергии связи водорода, так и расстояния от молекулы H2 до монослоя. Если указанную процедуру моделирования в АО-базисе выполнить не до конца, опустив, например, поправки Бойса–Бернарди, тогда по сравнению с ПВ-расчетами энергия связи будет переоценена почти в 2 раза, а расстояние – недооценено примерно на 5%. Хотя такие изменения в геометрии незначительно влияют на энергию связи, излишняя “компактность” водородных комплексов может привести к некорректным оценкам водородной емкости материала. Полученный оптимизированный базис для монослоя TOC может быть использован далее для исследования других свойств этого материала, а также свойств возможных его модификаций, возникающих, например, путем декорирования поверхности слоя атомами сторонних элементов или нанесения структурных дефектов.
Различные поправки к функционалу электронной плотности на дисперсионное взаимодействие, упомянутые выше, дают близкие оценки для энергии связи. С осторожностью стоит использовать только поправки Гримме третьего поколения при моделировании сорбции на чистых углеродосодержащих наноматериалах, т.к. можно получить завышенные оценки Ebind для некоторых конфигураций водорода. Таким образом, нам удалось согласовать предсказания энергии и длины связи молекулы H2 с монослоем ТОС, сделанные с помощью АО- и ПВ-базисов, однако оценка энергии связи оказалась неутешительной (60–90 мэВ/молекулу H2) и далекой от границы коммерческой привлекательности материала для применения в водородной энергетике. Предстоит выяснить, смогут ли это сделать модификации тетраоксо[8]циркулена.
БЛАГОДАРНОСТИ
Работа выполнена при финансовой поддержке Министерства науки и высшего образования Российской Федерации (государственное задание № FENU-2023-0011 Фундаментальные основы безопасных водородных технологий).
Конфликт интересов. Авторы заявляют, что у них нет конфликта интересов.
About the authors
E. V. Anikina
Institute of Natural Sciences and Mathematics, South Ural State University
Author for correspondence.
Email: anikate@inbox.ru
Russian Federation, 454080, Chelyabinsk
D. V. Babailova
Institute of Natural Sciences and Mathematics, South Ural State University
Email: anikate@inbox.ru
Russian Federation, 454080, Chelyabinsk
M. S. Zhilin
Institute of Natural Sciences and Mathematics, South Ural State University
Email: anikate@inbox.ru
Russian Federation, 454080, Chelyabinsk
V. P. Beskachko
Institute of Natural Sciences and Mathematics, South Ural State University
Email: anikate@inbox.ru
Russian Federation, 454080, Chelyabinsk
References
- Simanullang M., Prost L. // Int. J. Hydrogen Energy. 2022. V. 47. № 69. P. 29808. https://www.doi.org/10.1016/j.ijhydene.2022.06.301
- Wang Y., Yang P., Zheng L., Shi X., Zheng H. // Energy Storage Materials. 2020. V. 26. P. 349. https://doi.org/10.1016/j.ensm.2019.11.006
- Anglada E., Soler J.M., Junquera J., Artacho E. // Phys. Rev. B. 2002. V. 66. № 20. P. 205101. https://www.doi.org/10.1103/PhysRevB.66.205101
- Ferre-Vilaplana A. // J. Chem. Phys. 2005. V. 122. № 10. P. 104709. https://www.doi.org/10.1063/1.1859278
- Vilela Oliveira D., Laun J., Peintinger M.F., Bre- dow T. // J. Comput. Chem. 2019. V. 40. № 27. P. 2364. https://www.doi.org/10.1002/jcc.26013
- Begunovich L.V., Kuklin A.V., Baryshnikov G.V., Valiev R.R., Ågren H. // Nanoscale. 2021. V. 13. № 9. P. 4799. https://www.doi.org/10.1039/D0NR08554E
- Fritz P.W., Chen T., Ashirov T., Nguyen A.D., Dincă M., Coskun A. // Angew. Chemie. 2022. V. 61. № 17. P. e202116527. https://www.doi.org/10.1002/anie.202116527
- Baryshnikov G.V., Minaev B.F., Karaush N.N., Minae-va V.A. // RSC Adv. 2014. V. 4. № 49. P. 25843. https://www.doi.org/10.1039/c4ra02693d
- Karaush-Karmazin N., Baryshnikov G., Minaeva V., Panchenko O., Minaev B. // Comput. Theor. Chem. 2022. V. 1217. P. 113900. https://www.doi.org/10.1016/j.comptc.2022.113900
- Artacho E., Anglada E., Diéguez O., Gale J.D., García A., Junquera J., Martin R.M., Ordejón P., Pruneda J.M., Sánchez-Portal D., Soler J.M. // J. Phys. Condens. Matter. 2008. V. 20. № 6. P. 064208. https://www.doi.org/10.1088/0953-8984/20/6/064208
- Kresse G., Joubert D. // Phys. Rev. B. 1999. V. 59. № 3. P. 1758. https://www.doi.org/10.1103/PhysRevB.59.1758
- Soler J.M., Artacho E., Gale J.D., García A., Junque- ra J., Ordejón P., Sánchez-Portal D. // J. Phys. Condens. Matter. 2002. V. 14. № 11. P. 2745. https://www.doi.org/10.1088/0953-8984/14/11/302
- Salahdin O.D., Sayadi H., Solanki R., Parra R.M.R., Al-Thamir M., Jalil A.T., Izzat S.E., Hammid A.T., Arenas L.A.B., Kianfar E. // Appl. Phys. A. 2022. V. 128. № 8. P. 703. https://www.doi.org/10.1007/s00339-022-05789-2
- Grimme S. // J. Comput. Chem. 2006. V. 27. № 15. P. 1787. https://www.doi.org/10.1002/jcc.20495
- Grimme S., Antony J., Ehrlich S., Krieg H. // J. Chem. Phys. 2010. V. 132. № 15. P. 154104. https://www.doi.org/10.1063/1.3382344
- Perdew J.P., Burke K., Ernzerhof M. // Phys. Rev. Lett. 1996. V. 77. № 18. P. 3865. https://www.doi.org/10.1103/PhysRevLett.77.3865
- Dion M., Rydberg H., Schröder E., Langreth D.C., Lundqvist B.I. // Phys. Rev. Lett. 2004. V. 92. № 24. P. 246401. https://www.doi.org/10.1103/PhysRevLett.92.246401
- Berland K., Hyldgaard P. // Phys. Rev. B. 2014. V. 89. № 3. P. 035412. https://www.doi.org/10.1103/PhysRevB.89.035412
- Abinit’s pseudo database (2023) Fritz-Haber-Institute (FHI). https://departments.icmab.es/leem/SIESTA_MATERIAL/Databases/Pseudopotentials/periodictable-intro.html. Cited 20 May 2023.
- Созыкин С.А., Бескачко В.П., Вяткин Г.П. // Вестник ЮУрГУ. Серия “Математика. Механика. Физика”. 2015. T. 7. № 3. C. 78.
- Boys S.F., Bernardi F. // Mol. Phys. 1970. V. 19. № 4. P. 553. https://www.doi.org/10.1080/00268977000101561
- Anikina E.V., Beskachko V.P. // Bull. South Ural State Univ. Ser. “Mathematics. Mech. Physics.” 2020. V. 12. № 1. P. 55. https://www.doi.org/10.14529/mmph200107
- Artacho E., Sánchez-Portal D., Ordejón P., García A., Soler J.M. // Physica Status Solidi B. 1999. V. 215. Iss. 1. P. 809. https://www.doi.org/10.1002/(SICI)1521-3951 (199909)215:1<809::AID-PSSB809>3.0.CO;2-0
- Anikina E., Naqvi S.R., Bae H., Lee H., Luo W., Ahuja R., Hussain T. // Int. J. Hydrogen Energy. 2022. V. 47. № 19. P. 10654. https://www.doi.org/10.1016/j.ijhydene.2022.01.126
Supplementary files
