Ab initio molecular dynamics simulation of the superionic state in Pb_0.78Sr_0.19K_0.03F_1.97 solid solution: fluoride sublattice behaviour

Full Text

ВВЕДЕНИЕ

Интерес к фторидным материалам типа MF₂ (M = Ca, Sr, Ba) и гетеро- и гомовалентным твердым растворам на их основе обусловлен как их возможным использованием в электрохимических устройствах (источники тока, сенсоры, фтор-ионные насосы, оптоионные приборы) [1–3], так и применением их в качестве оптически активных материалов [4], перспективных функциональных покрытий [5]. Одним из модельных материалов для исследования особенностей фтор-ионного переноса является фторид свинца и твердые растворы на его основе [6–10].

Компьютерный эксперимент, в частности метод молекулярной динамики (МД) [11], является современным мощным инструментом, позволяющим изучать и прогнозировать характеристическое поведение кристаллов на атомном уровне. В классической МД для успешного воспроизведения экспериментальных зависимостей необходимо подобрать вид и параметры потенциала межатомного взаимодействия. Выбор параметров силового поля осуществляется на основе различных данных исследуемых объектов (структурные, упругие, термодинамические, диэлектрические и др.), что позволяет анализировать природу заранее выбранных характеристик, например подвижность ионных носителей в суперионных материалах. В то же время при использовании квантово-химических подходов (ab initio молекулярная динамика или неэмпирическая молекулярная динамика – НЭМД) взаимодействие между атомами определяется расчетами электронной структуры системы на каждом шаге движения атомов, т.е. без привязки к экспериментальным данным.

Известен ряд работ по моделированию фторида свинца с использованием классических МД-расчетов, основанных на модели жестких ионов. В [12–14] для PbF₂ был описан так называемый фарадеевский (“размытый”) фазовый переход в суперионное состояние (при температуре ~700 K). Экспериментально наблюдаемую температурную зависимость коэффициента диффузии в объемном кристалле удалось воспроизвести [15] при использовании большой расчетной ячейки из 1500 атомов и показать, что резкое возрастание проводимости обусловлено увеличением концентрации подвижных анионов. Фазовый α→β-переход из орторомбической в кубическую фазу, наблюдаемый при наложении давления, был изучен в [16].

Для адекватного воспроизведения термодинамических и тепловых свойств кристалла фторида свинца, как было показано в [17–19], необходимо учитывать влияние температуры на параметры межатомных потенциалов от температуры. Для корректного описания суперионного перехода температурный диапазон от 300 до 1000 К пришлось разбить на три участка: рассматривались две температурные области до температуры фарадеевского фазового перехода при T_фп = 700 K и одна – выше Т_фп. Можно утверждать, что результаты расчетов МД-методом позволили проанализировать разупорядоченное состояние фторида свинца (выше Т_фп) и показать важность наличия дефектов для движения атомов фтора в суперионном состоянии [13, 14, 19, 20].

Процессы плавления, кристаллизации и стеклообразования [20–22] в PbF₂ также исследовались с использованием метода МД. Было показано, что в расплаве и стекле наиболее вероятным координационным числом для атомов свинца является 8, что соответствует движению атомов фтора по локальным искажениям флюоритовой кристаллической решетки.

Состояние с высокой ионной проводимостью для твердых растворов на основе PbF₂ анализировалось в [23, 24], была подчеркнута определяющая роль допантов в повышении подвижности анионов фтора.

Таким образом, применяя метод классической МД, удалось успешно описать поведение кристалла фторида свинца в суперионном состоянии. Поэтому логично ожидать использование ab initio МД для анализа суперионного состояния фторидных материалов, тем более что для ряда простых фторидов MF₂ и LaF₃ ранее были осуществлены [25–28] расчеты электронной структуры. Отметим, что подход НЭМД эффективно применялся для изучения транспортных свойств в кислород-ионных проводниках на основе диоксида циркония [29–32] и ряда катионпроводящих твердых электролитов [33–36].

В настоящей работе для исследования динамических и структурных свойств и особенностей подвижности анионов фтора в твердых растворах PbF₂–SrF₂–KF использовали метод ab initio МД.

Выбор твердого раствора с двойным гомо- и гетеровалентным замещением атомов свинца обусловлен наличием экспериментальных результатов, которые указывают на высокие значения ионной проводимости в такой системе при комнатной температуре [7]. Введение в низкосимметричную матрицу a-PbF₂ небольшого количества (20 мол. %) SrF₂ с относительно высокой температурой плавления (1473°С) стабилизирует при комнатной температуре флюоритовую структуру образующегося твердого раствора. В свою очередь стабилизация фазы приводит к увеличению транспорта фтора, осуществляемому по междоузельному механизму, хорошо изученному для кристаллов со структурой флюорита. Дополнительное добавление гетеровалентного атома калия (5 мол. %) вместо атома свинца создает вакансию атомов фтора в анионной подрешетке, что еще больше повышает их подвижность. Полученный твердый раствор, отвечающий составу Pb_0.75Sr_0.20K_0.05F_1.95, обладает наивысшей ионной проводимостью в этой системе [7].

КОМПЬЮТЕРНЫЙ ЭКСПЕРИМЕНТ

Построение расчетной ячейки. Базовый кристалл PbF₂ имеет (в суперионном состоянии) достаточно простую кубическую структуру типа флюорита, пр. гр. Fm3̅m, с параметром элементарной ячейки a = 5.973 Å. Для МД-моделирования был создан расчетный бокс размером 2 × 2 × 2 элементарных ячеек вдоль кристаллографических направлений a, b и с соответственно, т.е. расчетный бокс состава Pb₃₂F₆₄ имел длину ребра куба 11.855 Å и объем 1666 ų. В настоящей работе помимо базового состава PbF₂ изучен твердый раствор состава Pb_0.78Sr_0.18K_0.03F_1.96, образованный случайным замещением атомов свинца атомами стронция и калия, что соответствует расчетной ячейке Pb₂₅Sr₆KF₆₃. Расчетный бокс для указанной системы показан на рис. 1.

 

Рис. 1. Расположение атомов в моделируемой системе Pb₂₅Sr₆KF₆₃ в начале (а) и после окончания (б) НЭМД-расчетов

 

МД-расчеты. Для проведения ab initio МД-расчетов использовали пакет программ Materials Studio и модуль DMol³, в котором реализован метод DFT на атомном базисе. Применяли базис DND, который соответствует базисному набору Gaussian 6-31G*, и обобщенный градиентный функционал PBE как один из “классических” функционалов – такое сочетание представляет собой приемлемый компромисс между точностью и необходимыми значительными компьютерными ресурсами [37–39].

НЭМД-расчеты проводили при модельной температуре 1000 K, поскольку в этом случае катионы металлов остаются малоподвижными и “держат” кристаллический остов, тогда как анионы фтора обладают значительными длинами свободного пробега в кристалле, что соответствует суперионному состоянию. Шаг интегрирования уравнений движения по времени был выбран 1 × 10⁻¹⁵ с (1фс), что обеспечивало стабильность полной энергии системы с точностью до 0.5 эВ. Длительность расчетов составляла 5 пс.

Подвижность частиц в кристаллической решетке характеризовали с помощью функций среднеквадратичных смещений (СКС):

$r_k^2\left(t\right)=\frac{1}{N}{\sum }_{k=1}^N\left(\begin{array}{c}{\left[x_k\left(t\right)-x_k\left(0\right)\right]}^2+\\ +{\left[y_k\left(t\right)-y_k\left(0\right)\right]}^2+{\left[z_k\left(t\right)-z_k\left(0\right)\right]}^2\end{array}\right)$, (1)

где x_k(t), y_k(t) и z_k(t) – координаты частиц сорта k в момент времени t. Временные зависимости  (за исключением начального участка) аппроксимировали прямыми, из величины тангенса угла наклона которых рассчитывали коэффициенты диффузии ионов (D_k), используя соотношение Эйнштейна:

$r_k^2\left(t\right)=6D_{k}t+B_k$, (2)

где коэффициент B_k описывает тепловые колебания частиц около положения равновесия.

Для анализа структурных особенностей фтор-ионной подсистемы изучали радиальные парные корреляционные функции (РПКФ):

$g_{ij}\left(r\right)=\frac{n_{ij}\left(r\right)}{\rho \cdot 4\pi r^{2}dr}$, (3)

где n – среднее число частиц j-го типа в сферическом слое толщиной dr на расстоянии r от частицы типа i;  – среднее значение атомной плотности моделируемой системы.

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

Рассчитанные РПКР распределения для системы Pb₂₅Sr₆KF₆₃ представлены на рис. 2. Как видно из рис. 2а, максимумы на кривой для пары Pb–Pb отвечают максимумам для идеальной флюоритовой решетки, в то время как максимумы на кривых для пар F–F и Pb–F существенным образом смещены относительно максимумов для идеального кристалла.

 

Рис. 2. Радиальные парные функции распределения атомов: пары F–F, Pb–Pb, Pb–F в системе Pb₂₅Sr₆KF₆₃ (а). Вертикальные прямые отвечают ситуации идеального кристалла с флюоритовой структурой; пары Pb–F (1) в кристалле β-PbF₂, пары Pb–F (2), Sr–F (3), K–F (4) в системе Pb₂₅Sr₆KF₆₃ (б)

 

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

На рис. 2б детально показаны РПКФ для всех пар катион–анион для первой координационной сферы. Видно, что для разных катионов расположение максимумов существенно разное. Для объяснения этого эффекта были рассчитаны усредненные заряды катионов, результаты расчетов приведены в табл. 1.

 

Таблица 1. Усредненные заряды на атомах (по Малликену) до и после выполнения расчетов и их ионные радиусы

Катион	Заряд до начала расчета, e	Заряд по окончании расчета, e	Ионный радиус, Å
Pb	1.202	1.150	0.119
Sr	1.800	1.782	0.118
K	0.911	0.968	0.138

 

Атом калия, имеющий минимальный абсолютный заряд (табл. 1) и наибольший ионный радиус по сравнению с атомами стронция и свинца, дает возможность атомам фтора максимально отойти от атома металла, что приводит к существенному смещению пика на РПКФ вправо относительно реперной точки – пика Pb–F в кристалле PbF₂ и уширению пика. Положение максимума РПКФ для пары Sr–F отвечает положению первой координационной сферы для идеального кристалла SrF₂; это означает, что атомы стронция пытаются организовать вокруг себя подобие самостоятельной фазы “нанокристаллов” SrF₂.

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

Поскольку моделируемая ячейка имеет сравнительно небольшие размеры, был проведен анализ локальной диффузии атомов фтора вблизи атомов металлов, т.е. в окрестности первой координационной сферы. Оказалось, что в процессе моделирования равновесие в системе достигается в течение ~1 пс, о чем свидетельствуют временные данные СКС, представленные на рис. 3. Действительно, на графике зависимости СКС от времени угол наклона апроксимирующей прямой практически не изменяется с увеличением времени, что соответствует установившемуся равновесному состоянию исследуемой системы.

 

Рис. 3. Среднеквадратичные смещения атомов фтора, локализованные в первой координационной сфере атомов металлов: Pb (1) в β-PbF₂, Pb (2), Sr (3), K (4) в системе PbF₂–SrF₂–KF

 

Коэффициенты самодиффузии атомов фтора были рассчитаны по данным среднеквадратичных смещений по формуле Эйнштейна (2) и приведены в табл. 2.

 

Таблица 2. Коэффициенты самодиффузии атомов фтора в кристалле β-PbF₂ и твердом растворе PbF₂–SrF₂–KF

Система	DF, см2/с
β-PbF2 (F около атомов Pb)	2.91 × 10−5
PbF2–SrF2–KF (F около атомов Pb)	4.00 × 10−5
PbF2–SrF2–KF (F около атомов Sr)	2.32 × 10−5
PbF2–SrF2–KF (F около атомов K)	3.19 × 10−5

 

Как видно из результатов расчетов, наименьшую подвижность имеют атомы фтора, локализованные около атомов Sr. Как отмечалось выше, катионы Sr пытаются организовать около себя нанокристаллы SrF₂, а подвижность фтора в образцах SrF₂ существенно ниже, чем в PbF₂ [6]. Это связано с тем, что высокая степень ионности Sr (табл. 1) обусловливает значительное кулоновское притяжение атомов фтора и тем самым затрудняет их транспорт. Значительная подвижность атомов фтора вблизи катионов калия может быть связана как со слабым зарядом К⁺ (табл. 1) и, соответственно, меньшей степенью взаимодействия с ионами фтора, так и наличием дополнительных вакансий как локального компенсатора недостающего заряда при замещении свинца на калий Pb→K. Однако наибольшая подвижность атомов фтора наблюдается в окрестности атомов свинца, что отвечает наблюдаемым экспериментальным данным и, видимо, связано с наиболее высокой электронной поляризуемостью атомов свинца по сравнению с атомами стронция и более сильной ионной связью между атомами стронция и фтора [7].

ЗАКЛЮЧЕНИЕ

Таким образом, метод НЭМД позволяет оценить подвижность атомов фтора в многокомпонентном твердом растворе с помощью расчетов электронной структуры, без привлечения эмпирических данных, что необходимо в случае метода классической МД. При этом данные расчетов методом НЭМД адекватно объясняют экспериментальные данные. Это дает возможность детального анализа как динамических процессов, происходящих при взаимодействии различных гомо- и гетеровалентных допантов и атомов фтора в базовой матрице исходного твердого электролита, так и и структурных особенностей, характерных для сложных нестехиметрических систем на атомном уровне.

Работа выполнена при использовании программно-аппаратных ресурсов ресурсного центра “Вычислительный центр” Санкт-Петербургского государственного университета, а также при поддержке Российского научного фонда (гранта № 22-23-00465). А.К. Иванов-Шиц выполнял работу по Госзаданию ФНИЦ «Курчатовский институт».

About the authors

A. V. Petrov

St. Petersburg State University

Author for correspondence.
Email: a.petrov@spbu.ru
Russian Federation, St. Petersburg

Q. Ji

St. Petersburg State University

Email: a.petrov@spbu.ru
Russian Federation, St. Petersburg

I. V. Murin

St. Petersburg State University

Email: a.petrov@spbu.ru
Russian Federation, St. Petersburg

A. K. Ivanov-Schitz

Shubnikov Institute of Crystallography of Kurchatov Complex of Crystallography and Photonics of NRC “Kurchatov Institute”

Email: a.petrov@spbu.ru
Russian Federation, Moscow

References

Supplementary files