Содержание к диссертации
Введение
Глава 1. Обзор литературы 11
1.1. Кристаллическая структура двойных фторидов редкоземельных элементов 11
1.2. Экспериментальные методы исследования двойных фторидов редкоземельных элементов под давлением 12
1.3. Исследование двойных фторидов редкоземельных элементов под давлением с помощью моделирования из первых принципов 28
1.4. Теория фазовых переходов Ландау 32
1.5. Заключение 36
Глава 2. Квантовомеханические методы описания твердого тела и их реализация 38
2.1. Теория функционала плотности 38
2.2. Заключение 56
Глава 3. Теоретические методы расчета механических и колеба тельных свойств твердых тел 57
3.1. Теория упругости 57
3.2. Уравнения состояния. Уравнение Берча-Мурнагана 66
3.3. Фононная дисперсия 69
3.4. Заключение 71
Глава 4. Результаты моделирования свойств соединений GdLiF4 и LuLiF4 73
4.1. Методика вычислений 73
4.2. Структурные свойства 78
4.3. Механические свойства 82
4.4. Электронные свойства 85
4.5. Колебательные свойства 89
4.6. Фазовые переходы 94
4.7. Заключение 103
Заключение 105
Список работ, опубликованных по теме диссертации 107
Литература
- Экспериментальные методы исследования двойных фторидов редкоземельных элементов под давлением
- Теория функционала плотности
- Уравнения состояния. Уравнение Берча-Мурнагана
- Механические свойства
Экспериментальные методы исследования двойных фторидов редкоземельных элементов под давлением
Давления, достигаемые в ячейке с алмазными наковальнями, могут привести к уменьшению объема более, чем в два раза, что, в свою очередь, может вызвать перестройку взаимного расположения атомов решетки с одновременным изменением межатомных расстояний [13]. При этом изменяется число соседних атомов, окружающих любой атом решетки (изменение координационного числа). Твердое тело с новым расположением атомов (с новой структурой) обладает и новыми, зачастую кардинально отличающимися свойствами по сравнению с исходной структурой. Таким образом, с помощью ячейки с алмазными наковальнями можно исследовать свойства фаз высокого давления, которые могут значительно отличаться от свойств фаз того же вещества при нормальных условиях.
Многочисленные современные экспериментальные исследования, в том числе Грзечника и др. [5, 8, 14–20], показывают, что оксошеелиты, чаще всего, переходят из тетрагональной структуры шеелита в моноклинную структуру вольфрамита (P2/c, Z=2) вблизи давления 11 ГПа; флуорошеелиты, с другой стороны, трансформируются несколько иначе – из той же структуры в моноклинную модификацию фергусонита (C2/c, Z=4 выше 300 K и P21/c, Z=4 от 0 до 300 К). Однако, обратимость в первоначальную структуру после малого гистерезиса наблюдается в шеелитах обоих видов [21].
Подробнее остановимся на флуорошеелитах. Грзечник и др. изучили структуры LuLiF4 [5], YLiF4 [17] и GdLiF4 [8] методом рентгеноструктурного анализа в ячейке с алмазными наковальнями при комнатной температуре под давлением. Модификация трёх структур происходила вблизи давления 11 ГПа, но с различными последствиями. В случае LuLiF4 наблюдался обратимый фазовый переход в структуру фергусонита (С12/с1, Z=4) при 10.7 ГПа, в то время как для YLiF4 переход при 10 ГПа приводил к структуре типа фергусонит c другой симметрией (I2/a, Z=4). GdLiF4 разлагался на твёрдые составляющие Gd1-LiF3-2 (P63/mmc, Z = 2) и LiF при давлении 13.1 ГПа.
Из-за различного поведения данных соединений под давлением стоит подробнее рассмотреть основные механизмы фазовых переходов в них. Соединение YLiF4 изоморфно соединениям с химической формулой LiF4 (=Eu-Lu) и имеет структуру шеелита [22, 23]. Данный кристалл был широко изучен при различных температурных условиях для лучшего понимания химических связей и взаимодействий между центральными катионами и их кристаллическим окружением [24, 25]. Что касается структурного поведения при высоких давлениях, измерение комбинационного рассеяния при комнатной температуре указывает на структурный фазовый переход при давлении 7 ГПа [26]. Хотя эксперименты по люминесценции на YLiF4:Eu при T=300 K не показывают какой-либо аномалии при 7 ГПа, обратимый фазовый переход ожидается в диапазоне давлений от 10 до 10.5 ГПа, и связан он с понижением симметрии кристалла [27]. Низкотемпературные измерения люминесценции в YLiF4, допирован-ном Nd, указывают на очень слабую аномалию структурных параметров фазы шеелита вблизи 6 ГПа и, кроме того, явно – на структурный переход вблизи 10.3 ГПа [28].
Непосредственно, в исследованиях Грзечника и др., был использован кристалл YLiF4, допированный 2 % Nd, как и в работе [28] по изучению люми 15 несценции. Дифракционные диаграммы YLiF4 были получены при различных давлениях при температуре T=300 K (см. Рис. 1.3). При давлениях до 10.5 ГПа все дифракционные пики соответствуют структуре шеелита.
Дифрактограммы порошка YLiF4 при различном давлении без предварительного отжига образца. Фон, возникающий в основном из-за комптоновского рассеяния в алмазной ячейке, вычитается из всех диаграмм [17].
При высоких давлениях разделение пиков и появление дополнительных слабых откликов указывают на фазовый переход в полиморфную модификацию со структурой, которая является неупорядоченной при низком давлении. Данная фаза высокого давления плавно переходит в другую полиморфную модификацию, начиная с давления 17 ГПа. Первая и вторая фазы высокого давления сосуществуют в широком диапазоне температур. При снятии давления, начиная от 24.8 ГПа, характерные черты дифракционных картин второй фазы наблюдаются вплоть до 10 ГПа. Фаза шеелита образуется вновь при полностью снятом давлении.
Экспериментальные исследования подтверждают, что первая фаза высокого давления имеет тип фергусонит (I2/a, Z=4), как и та, которая была впервые найдена в ортогалогенидах. В частности, данный результат был получен также и для редкоземельных ортониобатов и ортотанталатов [29–31]. Особенности ферроэластичных фазовых переходов под влиянием температуры в этих соединениях схожи с особенностями перехода типа шеелит-фергусонит в YLiF4 под давлением. В течение фазового перехода шеелит-фергусонит объем ячейки, а также тетрагональный и моноклинный параметры решетки изменяются плавно. Это указывает на то, что разница объемов между двумя полиморфными модификациями либо незначительна, либо ниже возможного предела измерения. Разность аксиальных соотношений / и / является мерой моноклинного разупорядочения и увеличивается с давлением.
Зависимости ( ) для полиморфных шеелита и фергусонита были аппроксимированы уравнением состояния Берча-Мурнагана вплоть до давления 10.7 ГПа. Результат аппроксимации представлен в сравнительной таблице 1.1, которая содержит экспериментальные значения механических параметров (модуль сжатия, первая производная модуля сжатия по давлению и равновесный объем) соединений YLiF4, GdLiF4, LuLiF4.
Теория функционала плотности
Пространственное расположение R& (совместно с соответствующими граничными условиями) отличает по существу молекулу от твердого тела. Подобным образом, слагаемое U отличает одно-частичную квантовую механику (см. уравнение (2.1) от крайне сложной многочастичной задачи, поставленной уравнением (2.2). Эти свойства заложены в методе DFT фундаментальным образом.
Обычный квантово-механический подход решения уравнения Шредингера (УШ) может быть записан как следующая последовательность: v(x) = гр(гі,Г2, ...,гдг) = наблюдаемые, (2.7) т.е. сначала система определяется посредством выбора (г), происходит подстановка потенциала в УШ, решается уравнение для волновой функции ф и затем вычисляются наблюдаемые посредством взятия математических ожиданий операторов с этой волновой функцией. Одна из наблюдаемых, которая вычисляется подобным образом, это - плотность числа электронов:
В течение десятилетий работы над задачей многих тел было развито множество мощных методов для решения уравнения Шредингера. В физике, например, один из методов решения этой задачи – диаграммная теория возмущений (основанная на диаграммах Фейнмана [154] и функциях Грина [155]), в то время как в квантовой химии часто используется метод конфигурационного взаимодействия [152, 156], который основан на разложении многоэлектронной волновой функции по детерминантам Слейтера, где учитывается антисимметричность волновой функции относительно взаимных перестановок любых пар электронов. Также существует большое количество более специальных техник. Проблема, возникающая при использовании этих методов, состоит в том, что они требуют огромных вычислительных ресурсов: зачастую бывает просто не возможно применить их эффективно к большим и сложным системам. Еще не удавалось вычислить химические свойства молекулы, состоящей из 100 атомов только с помощью метода конфигурационного взаимодействия или электронную структуру реального полупроводника, используя только функции Грина. Метод DFT же обеспечивает жизнеспособную альтернативу, возможно менее точную, но намного более универсальную. Точность в методе DFT является относительным понятием. Воспроизводимость экспериментальных результатов методом DFT в конкретных приложениях зависит от качества приближенных функционалов плотности, задействованных при моделировании той или иной системы. Для малого числа частиц или системы, обладающей специальной симметрией, могут быть получены точные решения уравнения Шредингера, и ни один приближенный функционал не сможет конкурировать с точными решениями. Современные сложные функционалы плотности позволяют достичь высокой точности для реальных систем. Длина связи молекул может быть предсказана со средней ошибкой менее, чем 0.001 нм, параметры решетки твердых тел -со средней ошибкой менее, чем 0.005 нм и энергия молекул в пределах 0.2 эВ [157-160]. Метод DFT явно указывает, что нерелятивистские кулоновские системы отличаются только их потенциалом v(r) и дает рецепт для работы с универсальными операторами Т и U. Более того, метод DFT обеспечивает способ систематически отображать многочастичную задачу с U на одно-частичную без U. Все это сделано путем присвоения плотности числа частиц п(г), одной среди многих наблюдаемых, статуса ключевой переменной, на которой базируется вычисление других наблюдаемых. Данный подход формирует основу большинства расчетов электронной структуры в физике и химии. Большая часть из того, что мы знаем об электрических, магнитных и структурных свойствах материалов вычислено с помощью метода DFT, и вклад, который данный метод внес в на 44 уку о молекулах, был отражен в 1998 году присуждением Нобелевской премии по химии Вольтеру Кону [161] и отцу-основателю DFT Джону Поплу [162], который сыграл важную роль в реализации данного метода в вычислительной химии.
Подход, используемый в теории функционала плотности, может быть записан как последовательность: n(r)= (ri,..., Гдг)= /(г), (2.9) т.е., знание плотности числа частиц п(г) предполагает знание волновой функции и потенциала и, как следствие, всех других наблюдаемых. Хотя эта последовательность описывает концептуальную структуру метода DFT, она не отражает реальной картины того, что сделано с помощью нее в реальных приложениях, которые обычно развиваются по разному и не используют явно многочастичные волновые функции.
Основой теории функционала плотности являются две теоремы Хоэнберга-Кона [144] и вычислительная схема, предложенная Коном и Шэмом [163] несколько позже. В 1964 году Хоэнберг и Кон доказали, что полная энергия основного состояния является функционалом зарядовой плотности и этот функционал минимален для соответствующей плотности. Вариационный принцип метода DFT обеспечивает основу для изучения взаимодействующих электронных систем через понятие электронной плотности. Однако, для прямого применения теорем Хоэнберга-Кона требуется аппроксимация функционалов кинетической и обменно-корреляционной энергии, что приводит к теориям Томаса-Ферми, которые дают более слабое соответствие с экспериментальными данными [164]. Теоремы Хоэнберга-Кона, сами по себе, являются значительным концептуальным достижением, но они не приводят немедленно к практическому уcпеху в вычислениях.
Первая из них утверждает, что для основного состояния уравнение (2.6) может быть инвертировано: при данной плотности основного состояния щ(г) возможно вычислить соответствующую основному состоянию волновую функцию о(гь r2, rw). Это означает, что фо является функционалом щ. Следовательно, все наблюдаемые в основном состоянии есть тоже функционалы щ. Если фо может быть вычислена из щ и наоборот, обе функции эквивалентны и содержат в точности одинаковую информацию. С первого взгляда это кажется невозможным: как может функция одной векторной переменной г быть эквивалентна функции N векторных переменных гі,...,г/у. Как может одна произвольная переменная содержать ту же самую информацию как и N произвольных переменных. Принципиальный факт, который делает это возможным состоит в том, что знание щ(г) несет в себе скрытую информацию о большем, нежели чем о просто произвольной функции /(г).
Уравнения состояния. Уравнение Берча-Мурнагана
В последние годы особенно возрос интерес к исследованию упругости кристаллов под давлением [210]. Существуют два основных способа вычисления констант упругости с помощью моделирования материалов из первых принципов [211]. Первый подход основывается на анализе полных энергий деформированных состояний материала, а второй – на анализе изменений рассчитанных значений напряжений, возникающих при деформации. Данный подход “напряжение - деформация” был описан в 1983 году в статье Нильсена и Мартина [212], посвященной ab initio расчету напряжения и его применению к вычислению констант упругости второго, третьего и четвертого порядков. С тех пор появились дополнительные разработки и приложения в широком диапазоне физических, химических, материаловедческих, геофизических и минералогических журналов.
Основным законом упругой деформации является линейная зависимость между силой и вызываемой ею деформацией. Этот закон выполняется только для малых смещений и характеризует линейную область, в которой возвращающая сила линейна относительно смещений атомов. Закон Гука в тензорном виде для кристалла представляет собой линейную зависимость между компонентами тензора напряжений и компонентами тензора деформаций. Тензор упругости является тензором 4-го ранга и в общем случае определяется матрицей, состоящей из 99=81 компоненты. Для симметричных тензоров 2-го ранга производный тензор 4-го ранга имеет 66=36 компонент. Закон Гука в тензорном виде запишется в виде:
Он описывает связь компонентов j(i=1,6) малого напряжения, вызванного приложением малой деформации к релаксированной кристаллической структуре, с компонентами j( = 1,6). В этом уравнении у- - константы упругости релаксированной структуры. Здесь напряжение имеет смысл локально приложенного давления. После записи в матричной форме = С, уравнение (3.1) представляет собой линейную систему из шести уравнений с шестью переменными , включающую 21 неизвестную постоянную упругости Cij. Тензор упругости является симметричным и определяется в общем случае 21 независимой компонентой [213]. Число упругих констант сокращается в зависимости от симметрии кристалла и в тетрагональных кристаллах GdLiF4 и LuLiF4 уменьшается до 7.
Модуль MT интегрирован в интерфейс VASP. Работа модуля основывается на трех шагах. Во-первых, модуль MT анализирует симметрию ячейки и определяет требуемые направления деформации, чтобы в дальнейшем на основе этих данных можно было рассчитать константы упругости. Затем выполняется расчет энергии с помощью модуля VASP для каждой деформированной ячейки. И на последнем шаге модуль MT анализирует результаты и вычисляет константы упругости (константы упругости и матрицы упругих деформаций вместе с их собственными значениями и собственными векторами). Также можно вычислить модуль сжатия, скорость звука и термодинамические свойства (в рамках модели Дебая).
Более подробно опишем механизм, использующийся в модуле MT, для расчета констант упругости. В настоящее время, ah initio вычисление констант упругости при нормальном давлении не является проблемой. Однако, вопрос, касающийся способа, позволяющего рассчитать константы упругости при произвольном давлении, требует отдельного внимания. Рассмотрим кристалл, сжатый посредством изотропного давления P до плотности р\. При таком сжатии происходит малая (но не инфинитезимальная) однородная деформация кристалла, которая переводит каждую точку R неискаженной решетки Браве в новое положение R/ в деформированной решетке: Щ = / (% + eij)Rj- (3.2) з При однородной деформации параметры є - являются константами, не зависящими от R (cij = ji), где индексы i,j указывают декартовы координаты и каждый из них пробегает по трем значениям; Sij - символ Кронекера. Разложение внутренней энергии на единицу массы кристалла в соответствии с тензором деформации Лагранжа:
Механические свойства
Для сглаживания функции населенности электронных уровней была использована функция Гаусса с параметром размытия =0.1 эВ. Были рассмотрены различные методы, учитывающие влияние спина на структурные и механические свойства GdLiF4, а именно: немагнитный (магнитные моменты атомов не рассматриваются), спин-поляризованный и спин-орбитальный.
Остановимся подробнее на последних двух методах. Спин-поляризованный метод применяется для систем, которые содержат неспаренные электроны или молекулярные радикалы. В таких системах число электронов со спином вверх может отличаться от числа электронов со спином вниз. В начале 1970 годов фон Барт, Хедин и Гуннарссон [260], а также Лундквист [261] адаптировали обобщенную теорию функционала плотности для спин-поляризованных систем. Это привело к появлению теории функционала спиновой плотности в приближении локальной спиновой плотности (LSD – local spin density). В теории функционала локальной спиновой плотности (LSDF – local spin density functional) фундаментальными величинами являются электронная (r) и спиновая плотности (r). Спиновая плотность определяется как разность между плотностями электронов со спином вверх и со спином вниз:
В теории LSDF обменно-корреляционный потенциал электронов со спином вверх в общем случае отличается от потенциала электронов со спином вниз. Следовательно, эффективный потенциал в уравнениях Кона-Шэма (см. уравнение (2.13)) становится зависящим от спина. В атомных единицах Хартри уравнения Кона-Шэма в спин-поляризованной форме запишутся как: [— + V ff(e )]ifjf(r) = є ф (г), (7 =[ или (7 =-1, (4.3) V fJr) = VH + V ,[p(r), cr(r)], (4.4) где VH - это электростатическиий потенциал электронной зарядовой плотности и Vxc - обменно-корреляционный потенциал. Обменно-корреляционный потенциал в теории LSDF зависит от электронной и спиновой плотности, как видно из уравнения (4.4). Существуют два набора одночастичных волновых функций: один для электронов со спином вверх, а другой для электронов со спином вниз (каждый с соответствующими одноэлектронными собственными значениями). В случае, когда плотности электронов со спинами вверх и вниз равны, спиновая плотность зануляется, и теория LSDF становится идентичной подходу LDF.
Перейдем к методу, учитывающему спин-орбитальное взаимодействие. Электроны вблизи атомных ядер достигают столь большой кинетической энергии, что релятивистские эффекты становятся заметными, поэтому при вычислениях электронной структуры необходимо их учитывать. Возрастание релятивистской массы электронов вблизи ядра приводит к уменьшению плотности электронного заряда вокруг ядра по сравнению с нерелятивистским случаем. Подобные эффекты изменяют относительную энергию , , и -орбиталей, которая влияет на механизмы связи и энергии. Релятивистские эффекты приводят к спин-орбитальному разделению и возникает необходимость в использовании теории, учитывающей эти эффекты. Релятивистская теория электронной структуры основывается на решении уравнений Дирака и учитывает спин-орбитальное взаимодействие, что подробно обсуждается в книге Бьёркена, Дрелла [262] и Мессия.
Также в расчетах были использованы два способа рассмотрения f -электронов Gd: f - электроны заменены псевдопотенциалом (не включены в расчёт напрямую), и f - электроны учитываются в расчёте и рассматриваются как валентные. В дальнейшем, при исследовании структуры LuLiF4 магнитные моменты не рассматривались, и f -электроны Lu были включены в псевдопотенциал, так как методы, рассмотренные при изучении GdLiF4 (изоструктурное LuLiF4 соединение), не оказали значительного влияния на его структурные и механические свойства.
Далее была рассмотрена структура LuLiF4 с различными симметриями – I41/a и C12/c1 – в диапазоне давлений от 0 ГПа до 20 ГПа с шагом 2 ГПа. Следующие параметры были выбраны для получения энергетической сходимости: k-сетка содержала 9 9 9 точек в зоне Бриллюэна для структуры с симметрией I41/a и 7 7 9 точек для структуры с симметрией С12/с1. Для каждой k-точки волновые функции были разложены по базису плоских волн до максимальной энергии 796 эВ для симметрии I41/a и 576 эВ для симметрии С12/с1. Для сглаживания функции населенности электронных уровней при оптимизации сил и энергии применялась схема Метфесселя-Пакстона с параметром размытия 0.12 эВ для симметрии I41/a и 0.17 эВ для симметрии С12/с1. Данные параметры были выбраны таким образом, что дальнейшее их изменение не приводило к изменению рассчитываемых характеристик. Известно, что Lu имеет электронную конфигурацию 4f 145d16s2. Сильные корреляции, возникающие за счет учета d-орбиталей Lu, рассмотрены в приближении GGA+U [184]. Были выполнены расчеты для разных значений параметров Хаббарда в диапазоне от 1 до 3 эВ. Для дальнейших вычислений был подобран параметр = 1 эВ. Для расчета механических свойств был применён модуль MT. Параметр величины относительной механической деформации был равен 0.002. Колебательные свойства решетки были расчитаны в гармоническом приближении с помощью модуля PHONON [249]. Матрица силовых констант была получена на основе смещений атомов на величину ± 0.02 A. Энергия “отсечки” составила 500 эВ. В расчетах не учитывались продольный оптический (LO) и поперечный оптический (TO) режимы разделения мод.
Была рассмотрена структура GdLiF4 c начальными параметрами решетки = = 5.235(1) A, = 11.01(2) A, = 302 A3, которые были взяты из эксперимента [8]. Далее структура была оптимизирована в дипазоне давлений от 0 до 16 ГПа с шагом 2 ГПа, таким образом были найдены конфигурации с минимальной энергией. Изменение параметров решетки для разных методов в зависимости от давления представлена на рисунке 4.2 а), б). Также на данном рисунке представлены зависимости объема элементарной ячейки и аксиального отношения / от давления (Рис. 4.2 в), г)).
Графики показывают хорошее соответствие между экспериментальными данными [8] и расчетными значениями параметров (экспериментальные значения также представлены на рисунках). Основываясь на этих результатах, можно заключить, что учет спиновых степеней свободы не оказывает сильного влияния на структурные параметры для рассматриваемой системы. Наилучшее согласие с экспериментальными данными показывает параметр решетки в случае, когда -электроны включены в псевдопотенциал (см. Рис. 4.2 а)).