Содержание к диссертации
Введение
Глава 1. Расчет свойств электронной подсистемы в псевдопо тенциальном подходе 30
1.1. Исследование границ применимости псевдопотенциалов 30
1.2. Электронная теплоемкость 37
1.3. Тепловое давление электронной подсистемы 40
1.4. Расчет плотности состояний электронов из многочастичного расчета 41
1.5. Оценка степени ионизации плазмы ударно-сжатых металлов 46
1.6. Выводы к первой главе 49
Глава 2. Расчеты теплофизических свойств в квазигармониче ском приближении 52
2.1. Расчеты фононных свойств 52
2.2. Расчеты кривых плавления металлов с использованием критерия Линдемана 60
2.3. Влияние температуры электронной подсистемы на температуру плавления металлов 65
2.4. Выводы ко второй главе
Глава 3. Квантовые молекулярно-динамические расчеты ударно волновых экспериментов 73
3.1. Параметры квантового молекулярно-динамического расчета 73
3.2. Расчеты сплошных, пористых и повторных ударных адиабат Гюгонио 79
3.3. Реконструкция кривых изоэнтропического расширения 88
3.4. Выводы к третьей главе 100
Заключение 101
Приложение А. Давление и внутренняя энергия в алюминии, полученные из первопринципного квантового молекулярно динамического расчета 102
Литература
- Электронная теплоемкость
- Оценка степени ионизации плазмы ударно-сжатых металлов
- Расчеты кривых плавления металлов с использованием критерия Линдемана
- Реконструкция кривых изоэнтропического расширения
Электронная теплоемкость
Фундаментальная идея псевдопотенциальной теории — с помощью некоторой процедуры уменьшить степень осцилляции рассчитываемых валентных волновых функций вблизи ядра, тем самым уменьшив количество плоских волн, которые потребуются для их разложения. Более того, рассчитываются электроны только верхних оболочек, так как известно, что от поведения этих электронов зависит большинство физических свойств веществ. Орбита-ли внутренних электронов считаются неизменными — поведение волновых функций этих электронов не меняется при изменении внешнего химического окружения атома. Таким образом, суть данного подхода сводится к замене сильного электрон-ионного потенциала гораздо более слабым потенциалом -псевдопотенциалом - который бы описывал все основные свойства валентных электронов, включая релятивистские эффекты. Для этого исходная система заменяется системой, состоящей из псевдовалентных электронов и псевдоионного кора. Свойства псевдоиона таковы, что его потенциал вне некоторого радиуса обрезания (кора) совпадает с потенциалом истинного иона, но внутри этой сферы он гораздо слабее [13]. Тот факт, что потенциал внутри кора ослаблен — существенен, так как именно это делает решение уравнения Шре-дингера внутри сферы кора значительно проще. В этом случае осциляции ая
Схематичная иллюстрация полноэлектронного (сплошная линия) и псевдоэлектронного (штриховая линия) потенциалов и соответствующих им волноввіх функций. Радиус, на котором полноэлектронное и псевдопотенциалвное приближения совпадают, отмечен как гс волновых функций вблизи ядра уменьшаются или исчезают совсем. Ионный потенциал, валентная волновая функция и соответствующие псевдопотенциал и псевдоволновая функция изображены схематично на рис. 1.
В самом общем виде псевдопотенциал может быть представлен как где \lm) — сферические гармоники и V/ — псевдопотенциал углового момента /. Действуя на электронную волновую функцию, этот оператор раскладывает ее на сферические гармоники, каждая из которых умножается на соответствующий псевдопотенциал V/, индекс «NL» в данном случае означает «нелокальный». Большинство идей псевдопотенциального подхода пришло из метода ор-тогонализованных плоских волн (OPW), введенного Герингом [14, 15] в 1940 году, согласно которому базисные функции валентных состояний x qPW УД влетворяют уравнению: OPW / \
Метод ортогонализованных плоских волн стал основой для первых зонных расчетов веществ, отличных от sp-металлов [16]. Этот метод является предшественником современных псевдопотенциальных подходов и метода проекционных присоединенных волн (PAW).
Основные принципы, на которых базируется создание современных псевдопотенциалов, были введены в работе [17]. Согласно предложенной схеме, внутри некоторого атомного радиуса полноэлектронные волновые функции заменяются псевдоволновыми функциями, причем должен строго выполняться ряд условий: валентные псевдоволновые функции, генерируемые псевдопотенциалом, не должны содержать узлов; полноэлектронные и псевдовалентные собственные значения должны совпадать для выбранной в качестве прототипа атомной конфигурации; полноэлектронные и псевдоволновые функции должны совпадать вне выбранного радиуса кора гс; объемные интегралы по полноэлектронной и псевдозарядовой плотности внутри кора должны быть равны (правило сохранения нормировки); логарифмические производные полноэлектронной и псевдоволновой функции, а также их первые производные по энергии должны совпадать вне кора.
Псевдопотенциалы, удовлетворяющие этим условиям принято называть сохраняющими нормировку (norm-conserving).
Позже было установлено, что для хорошей переносимости псевдопотенциала, т.е. свойства валентных псевдоэлектронов описывать ответ системы на изменение химического окружения, требуется, чтобы кор располагался вокруг самого дальнего максимума полноэлектронной волновой функции, потому что только в этом случае распределение заряда и моменты полноэлектронной волновой функции хорошо воспроизводятся псевдоволновой функцией [18, 19]. Поэтому для элементов с сильно локализованными орбиталями (элементы первой группы, переходные металлы и редкоземельные элементы) получающийся псевдопотенциал требует большого набора базисных плоских волн. В качестве компромиссного решения часто идут на значительное увеличение радиуса кора, но такой подход неудовлетворительно сказывается на переносимости, и для любых новых химических сред требуются дополнительные тесты для проверки надежности таких мягких псевдопотенциалов.
Элегантное решение этой проблемы было предложено Вандербильтом [7]. Важной его особенностью является ослабление условия сохранения нормировки. Он предложил значительно сгладить волновые функции внутри кора и для компенсации дефицита заряда ввести локализованные в центре атома дополнительные заряды. Эти заряды обусловлены разницей между полноэлектронной и псевдоволновой функциями. Соответствующий псевдопотенциал из псевпотенциала сохраняющего нормировку полноэлектронная из ультра мягкого псевдопотенциала Радиальная волновая функция 2р орбитали кислорода. Сплошная линия — полноэлектронная волновая функция, точечная линия — псевдоволновая функция, сгенерированная с помощью псевдопотенциала, сохраняющего нормировку [17], штриховая линия — гладкая часть псевдоволновой функции, полученной с помощью ультрамягкого псевдопотенциала [7]. получил название ультрамягкий (US-PP). Радиус кора псевдопотенциала в таком случае может быть выбран как половина межчастичного расстояния, независимо от положения максимума полноэлектронной волновой функции, что продемонстрировано на рис. 2. Но для точного воспроизведения распределения зарядовой плотности и моментов полноэлектронной волновой функции требуется подгонка параметров вспомогательных зарядов.
Оценка степени ионизации плазмы ударно-сжатых металлов
Различные свойства металлов в кристаллической фазе необходимы как для решения фундаментальных задач, так и для огромного числа приложений. Современные вычислительные подходы физики твердого тела позволяют проводить весьма точные расчеты электронной структуры кристаллов и спектра колебаний решетки. В последнем случае традиционно применяется квазигармоническое приближение, в рамках которого, в частности, можно вычислить свободную энергию кристалла. Чаще всего этот метод применяется в упрощенной формулировке, в которой по фононному спектру при нулевой температуре электронов определяется температура Дебая, а далее используется выражение для свободной энергии в рамках классической модели Дебая.
Между тем, в последние годы этот метод получил активное развитие в связи с ростом вычислительных мощностей и появлением возможности последовательного учета всех мод колебаний кристалла. Пример такого расчета описан в статье [82].
Учет термического возбуждения электронов и теплового расширения кристалла существенно расширяет границы применимости данного подхода в области высоких температур вплоть до плавления. Благодаря тому факту, что знание свободной энергии определяет все термодинамические свойства системы, данный первопринципный метод может претендовать на построение полного термического уравнения состояния кристалла [83]. Современные исследования показывают, что при больших температурах вклад термического возбуждения электронов в электрон-фононное взаимодействие может становиться существенным, оказывая влияние на термодинамические свойства вещества. Подобное исследование приводится, например, в статье [84]. Особенно заметным вклад электронов становится в случае с нагретой электронной подсистемой, когда электронная температура существенно превосходит температуру кристалла. Экспериментальные исследования таких систем стали возможны лишь в последнее десятилетие с появлением мощных лазеров, способных производить ультракороткие импульсы.
В качестве исследуемых элементов были выбраны алюминий и медь. В расчетах использовался программный пакет VASP [23, 67-69] и код Phonopy [85, 86] для обработки первопринципных данных и получения фононных свойств.
Для первопринципных расчетов использовались PAW потенциалы [8, 24]. Для обменно-корреляционного функционала использовалось обобщенно-градиентное приближение (GGA) в параметризации Пердью, Бурке и Эрзенхофа (РВЕ) [71, 72]. Для расчета собственной энергии использовалось GW-прибли-жение [87]. Псевдопотенциал для алюминия содержал 3 валентных электрона, для меди — 11 электронов вне кора. Энергия обрезания плоских волн равнялась 700 эВ.
В ходе расчета силы, полученные из первопринципного подхода для супер-ячейки в прямом пространстве с использованием теории возмущений функционала плотности (DFPT)[88], реализованном в VASP, использовались для получения силовых констант и на их основе фононных частот с помощью кода Phonopy. Используемый в Phonopy метод [89] является модификацией подхода Парлински-Ли-Кавазу [90]. Для обоих металлов использовалась супер-ячейка ГЦК размером 2x2x2, содержащая 32 частицы, с решеткой k-точек в зоне Бриллюэна 8x8x8, генерируемой согласно схеме Монхор ста-Пака [73]. Для расчета фононных свойств использовалась сетка q-точек 61 X 61 X 61.
Свободная энергия кристаллической системы в приближении Борна-Оп-пенгеймера в рамках квазигармоничности представляется в следующем виде: F(V,T) = E0(V) + Fe(V,T) + Fph(V,T), (2.1) где EQ - энергия кристалла, Fe - свободная энергия электронной подсистемы, Fph - свободная энергия фононов.
Первые два слагаемых в уравнении (2.1) могут быть получены напрямую из первопринципного расчета. В свою очередь, фононный вклад в свободную энергию системы при постоянном объеме может быть получен из знания плотности состояния фононов, как функции фононных частот, посредством следующего выражения [83]:
Для получения равновесного значения объема кристалла при нормальных условиях проводилась серия расчетов для разных параметров решетки, после чего, используя выражение (2.1), исследовалась зависимость свободной энергии от объема для нахождения минимумов на изотермах. Результаты подобного исследования для алюминия продемонстрированы на рис. 2.1.
Для проверки используемого подхода и получаемых данных было проведено сравнение фононного спектра алюминия при нормальных условиях с экспериментальными данными. Результат сравнения представлен на рис. 2.2, где изображена кривая дисперсии фононов при 300 К.
На основе зависимостей свободной энергии кристаллов от температуры были получены кривые для коэффициента линейного расширения алюминия и меди, определяемого выражением (2.3), которые представлены на 14.5 64 66 68 70 72 74 76 Объем [AS ]
Рис. 2.1. Зависимость свободной энергии ГЦК-кристалла алюминия от объема вдоль изотерм. Синие линии с точками — кривые свободной энергии для разных температур от 0 до 1000 К (верхняя — 0 К, нижняя — 1000 К). Красная линия — минимумы энергии. Рис. 2.2. Кривая дисперсии фононов для алюминия при 300 К. Сплошная красная линия — данная работа, круглые маркеры — экспериментальные данные [91]. рис. 2.3-2.4. ад = зНЭ; (2 3)
Кроме того, были проведены дополнительные расчеты с другими псевдопотенциалами и обменно-корреляционными функционалами, чтобы проанализировать влияние данных факторов. Для алюминия проводилось две серии расчетов: с обобщенно-градиентным (GGA) и приближением локальной плотности (LDA) для обменно-корреляционного функционала; для меди — 3 серии: дополнительно к разным обменно-корреляционным функционалам были проведены расчеты для псевдопотенциала с большим числом валентных электронов вне кора. Сравнение проводилось с экспериментальными кривыми коэффициентов линейного расширения алюминия и меди, полученными усреднением имеющихся экспериментальных данных для выбранных металлов [92].
Из графиков 2.3-2.4 отчетливо видно, что структуры, описываемые в рамках обобщенно-градиентного приближения (GGA) обменно-корреляционного функционала, сильнее подвержены термическому расширению, нежели в приближении локальной плотности (LDA). В целом для исследуемых металлов до температуры 300 К лучше описывает экспериментальные данные обобщенно-градиентое приближение, хотя при дальнейшем увеличении температуры оно переоценивает расширение кристалла. В работе [93] проводилось комплексное первопринципное исследование коэффициентов теплового расширения для 4с! металлов, и был сделан вывод, что локальное приближение (LDA) обменно-корреляционного функционала дает более точные результаты до Т 300 К. В то же самое время, в работах [94, 95] было показано, что последовательный учет различных факторов позволяет улучшить совпадение результатов расчетов теплового расширения с экспериментом для обоих приближений.
Расчеты кривых плавления металлов с использованием критерия Линдемана
Квантовый молекулярно-динамический расчет существенно расширяет возможности исследования свойств плотной плазмы, ведь из такого расчета можно получить не только свойства электронной подсистемы, но и свойства ионов без привлечения дополнительных моделей.
Учитывая фундаментальную и прикладную важность изучения свойств ударно-сжатых веществ, в последней части работы была исследована возможность и реализованы методы по получению из квантовых молекулярно-дина-мических расчетов данных, сопоставимых по составу и качеству с результатами ударно-волновых экспериментов. В качестве исследуемого материала был выбран алюминий в виду наличия большого количества доступной для сравнения экспериментальной информации.
Квантовые молекулярно-динамические расчеты производились с использованием обобщенно-градиентного приближения (GGA) с коррекциями Пер-дью и Вонга (PW91) [5, 70] для обменно-корреляционного функционала, в качестве псевдопотенциала для алюминия был выбран ультрамягкий псевдопотенциал Вандербильта (US-PP) [7, 22]. Моделирование производилось для супер-ячейки, содержащей 108 атомов алюминия, с периодическими граничными условиями. Начальное расположение атомов соответствовало идеальной гранецентрированной (ГЦК) решетке. Плотность вещества р в расчете задавалась параметром решетки и оставалась постоянной. Моделирование динамики атомов алюминия производилось для временного интервала не менее 1 пс с шагом 2 фс. В расчетах для низких температур и плотностей время моделирования было увеличено до 8 пс, чтобы система смогла выйти на термодинамическое равновесие. Температура ионной подсистемы в расчете регулировалась посредством термостата Нозе-Гувера, температура электронной подсистемы задавалась через распределение Ферми-Дирака и равнялась ионной. Распределение электронной плотности в системе расчитывалось с помощью метода функционала плотности на каждом временном шаге. Затем по теореме Геллманна-Фейнмана рассчитывались силы, действующие на ионы, и определялись новые координаты и скорости атомов согласно уравнениям Ньютона. Для определенния зонной структуры использовалась одна точка (Г-точка) в зоне Бриллюэна. Энергия обрезания плоских волн равнялась 16 Ry. Термодинамические параметры системы определялись усреднением соответствующих величин после выхода системы на равновесие, статистическая погрешность результатов в расчетах была менее 1%. Кинетический вклад ионов в давление (pRT) добавлялся к величине (R - газовая постоянная), получаемой при моделировании.
Все расчеты были выполнены с псевдопотенциалом, содержащим 3 электрона вне кора, для диапазона температур 0.3-70 кК и плотностей 1-7.3 г/см3. Целесообразоность использования данного псевдопотенциала для рассмотренного диапазона температур и плотностей обсуждалась в Главе 1. При больших температурах начинается возбуждение электронов кора, что приводит к недооценке электронной теплопроводности [112], поэтому при Т 8 эВ необходимо использовать псевдопотенциал с большим числом валентных электронов.
Чтобы изучить влияние входных параметров на результаты моделирования был проведен ряд тестовых расчетов с большим значением энергии об резания плоских волн, большим числом атомов и /с-точек в зоне Бриллюэна, разными псевдопотенциалами и типом обменно-корреляционного функционала.
Как можно увидеть из рис. 3.1, выбор обменно-корреляционного функционала ведет к самым значительным различиям в результатах. Приближение локальной плотности (LDA) [113] для обменно-корреляционного функционала дает более низкие давление и энергию, нежели обобщенно-градиентное приближение (GGA). PAW потенциал с обменно-корреляционным функционалом (GGA) в параметризации Пердью, Бурке и Эрценхофа (РВЕ) [71, 72] дает до 3% большую энергию и меньшее давление, чем ультрамягкий псевдопотенциал (US-PP). Использование 256 частиц вместо 128, увеличение энергии обрезания плоских волн до 44 Ry, а также использование сетки 4x4x4 /с-точек в зоне Бриллюэна вместо Г-точки оказывает незначительное влияние на результаты расчетов.
Эволюция давления в системе в процессе расчета для разных временных шагов моделирования приведена на рис. 3.2. Из графика видно, что разогретая система быстро пришла в равновесие вне зависимости от выбранного временного шага.
Чтобы проверить качество выбранного ультрамягкого псевдопотенциала с 3-мя валентными электронами, были выполнены квантовые молекулярно-динамические расчеты изотермы 293 К.
Из графика 3.3 видно, что результаты расчета совпадают данными полноэлектронного расчета [114] и экспериментальными данными для ГЦК решетки [115, 116]. Таким образом, можно утверждать, что использование псевдопотенциала не вносит дополнительных ошибок в данном диапазоне давлений.
Реконструкция кривых изоэнтропического расширения
Сдвиг крайне правых экспериментальных точек на рис. 3.12 в область более высоких скоростей расширения может объясняться тем, что алюминий в данных точках находится в жидкой фазе [139], а, следовательно, может иметь место неоднородный разлет вещества в воздух. Для наглядности данного факта на рис. 3.15 приведены графики всех исследуемых изоэнтроп, рассчитанных согласно уравнению состояния [139], с указанием фазовых границ.
Еще две серии экспериментальных данных по изоэнтропическому расширению были получены при более высоких давлениях в ударно-сжатом веществе — 200-500 ГПа [123, 149]. В работе [123] в качестве среды с низким акустическим импедансом использовался воздух, а в работе [149] — кремниевый аэрогель.
На рис. 3.13 показаны изоэнтропы, соответствующие экспериментам [123]. Снова можно отметить хорошее согласие между результатами КМД и УРС расчетов и при этом сдвиг всех экспериментальных точек в область более высоких скоростей расширения. Нужно сказать, что эти экспериментальные точки не удается описать ни одним из существующих уравнений состояний. Анализ, проведенный в работе [139], показывает, что и начальные и конечные точки всех изоэнтроп [123] лежат в жидкой фазе при достаточно низкой плотности, поэтому такой результат может снова объясняться неоднородным разлетом поверхности вещества. Также обратим внимание, что полуэпириче-ское УРС предсказывает испарение алюминия при более низких плотностях, чем те, что были зафиксированы в экспериментах [123]. Вообще, равновесное испарение приводит к существенному росту скорости расширения вдоль изоэнтропы. Этот эффект заметен на результатах УРС на графике 3.13 для двух правых изоэнтроп, но давление испарения почти на порядок ниже, чем было получено в эксперименте для крайней правой изоэнтропы, а для других изоэнтроп эта разница еще больше. Более того, в работе [139] утверждается, 2 3 4
Изоэнтропы разгрузки в алюминии соответствующие экспериментам [123, 147— 149] в координатах давление-энтропия. Линии — расчеты с помощью полуэмпирического уравнения состояния [139]. На диаграмме отмечены фазовые состояния вещества. что алюминий не может испаряться при расширении в воздух после ударно-волнового сжатия (см. рис. 3.15).
В последней рассматриваемой экспериментальной работе [149] изучалось расширение алюминия в аэрогель. Результаты сравнения расчетов представлены на графике 3.14. Было получено отличное совпадение обеих расчетных моделей с экпериментальными измерениями. Разница давлений для начальных и конечных измеренных в эксперименте точек изоэнтроп менее одного порядка, а плотность жидкого алюминия при данных условиях даже выше, чем плотность твердого алюминия при нормальных условиях. Поэтому можно утверждать, что в этом эксперименте наблюдалась однородная поверхность расширяющегося в аэрогель алюминия [149], и автомодельное решение для адиабатического течения (3.4) может быть применено с хорошей точностью.
Интересно также было исследовать изменение температуры вдоль КМД-изоэнтроп в сравнении с результатами модели полуэмпирического УРС. На рис. 3.16 этот анализ продемонстрирован для нескольких изоэнтроп из экспериментов [123, 149]. Как было показано ранее, температуры на ударной адиабате, согласно обеим моделям, совпадают до области плавления, но при высоких плотностях УРС всегда выдает температуру выше, чем следует из КМД-расчета, и эта разница растет с увеличением плотности. Похожий эффект можно наблюдать и для изоэнтроп: широкодиапазонное УРС завышает температуру на изоэнтропах [123] по сравнению с КМД более чем на 20%. Это является следствием плохой калибровки полуэмпирической модели уравнения состояния в жидкой фазе по температурным данным.
В работе также была исследована изоэнтропическая скорость звука вдоль ударной адиабаты, которая измерялась в экспериментах [120, 150, 151]. Изоэнтропическая скорость звука позади фронта ударной волны (dP/dp)s может быть определена как производная вдоль изоэнтропы из точки на ударной
Скорость звука в алюминии за фронтом ударной волны. Красная кривая — данная работа, черная кривая — широкодиапазонное полуэмпирическое уравнение состояния [133], маркеры — экспериментальные точки [120, 150, 151]. адиабате в координатах давление-плотность. Для расчета производной вдоль изоэнтропы в заданной точке ударной адиабаты использовались полученные ранее интерполяционные функции E(V,T) и P(V,T). Результаты представлены на рис. 3.17. Наклон полученной кривой совпадает с графиком, полученным из УРС, а согласие с экспериментальными данными довольно хорошее. Немонотонное поведение кривой может объясняться погрешностью интерполяции в области плавления.
Результаты квантовых молекулярно-динамических расчетов для алюминия, использовавшиеся для интерполяции поверхностей энергии E(V, Т) и давления Р(у,Т): были сведены в таблицы и приведены в Приложении А.
Первопринципный метод КМД способен с высокой точностью описать все типы ударно-волновых экспериментальных данных: ударные адиабаты сплошного и пористого веществ, адиабату повторного сжатия, изоэнтропы разгрузки и скорость звука за фронтом ударной волны.
Сравнение с полуэмпирическим уравнением состояния показало отличное согласие результатов по давлению, плотности, скорости ударной волны и массовой скорости. Тем не менее, было отмечено различие в температурных данных для ударной адиабаты и изоэнтроп разрузки. Это может являться следствием недостатка надежных температурных измерений ударно-сжатого вещества для калибровки полуэмпирического уравнения состояния.
В целом КМД-расчет продемонстрировал хорошую точность расчета термодинамических свойств вещества, что дает основания для использования данного метода для калибровки полуэмпирических уравнений состояния в случае отсутствия экспериментальных данных.
