ПЛАВЛЕНИЕ И ДИНАМИКА РЕШЕТКИ НАТРИЯ ПРИ ВЫСОКИХ ДАВЛЕНИЯХ. РАСЧЕТ МЕТОДОМ КВАНТОВОЙ
МОЛЕКУЛЯРНОЙ ДИНАМИКИ
С. В. Лепешкина'ь* М. В. Магницкаяс, Н. Л. Мацко", Е. Г. Максимов а
а Физический институт им. П. Н. Лебедева Российской академии паук 119991, Москва, Россия
ьМосковский физико-технический институт Ц1700, Долгопрудный, Московская обл., Россия.
'"Институт физики высоких давлений Российской академии наук Ц2190, Троицк, Московская. обл., Россия.
Поступила в редакцию 8 ноября 2011 г.
Плавление и динамика решетки натрия исследованы при давлениях до 1 Мбар и температурах до 1000 К методом квантовой молекулярной динамики, т.е. с полным учетом энгармонизма. Результаты моделирования хорошо согласуются с экспериментом и с расчетом ab initio в квазигармоническом приближении, выполненным нами ранее. Расчеты показывают, что ангармонические взаимодействия оказывают малое влияние на кривую плавления и фононные частоты в Ка вплоть до температур вблизи плавления.
Е. Г. Максимов
1. ВВЕДЕНИЕ
В работе [1] в экспериментах на алмазных наковальнях при давлениях примерно до 140 ГПа было обнаружено необычное поведение кривой плавления Тт(р) натрия. Кривая Тт(р) проходит через ярко выраженный максимум в ОЦК-фазе (Ттаз, ~ 1000 К при р ~ 30 ГПа) и затем, в более плотной ГЦК-фазе резко падает до комнатных температур. Минимальное значение Г,„ , зафиксированное при р ~ 115 ГПа вблизи структурного перехода из ГЦК-фазы в менее симметричную структуру с 116, составляет около 300 К, что ниже точки плавления Na при нормальном давлении (примерно 371 К). Ранее аномальное плавление в столь широком интервале давлений и температур не наблюдалось ни в одном веществе и лишь недавно появилось сообщение о сходном поведении кривой плавления лития [2]. Опубликованные вскоре после работы [1] расчеты плавления Na [3 5], выполненные методом квантовой молекулярной динамики, с разной степенью точности воспроизводят экспериментальные данные. Например, в работе [3] получен максимум на кривой плавления и значительное уменьшение Тт в области перехода
*E-mail: lepeshkin'ffllpi.ru, lsw322'fflmail.ru
ГЦК —¥ с 116, при этом, как и в работе [4], основное внимание уделялось изучению локальной структуры жидкой фазы. В работе [5] исследовалось, в частности, смягчение упругих констант в твердой фазе, однако общее согласие полученных результатов с экспериментом оказалось неудовлетворительным, особенно при высоких давлениях. Следует отметить, что перечисленные работы не дают ясного представления о физических причинах столь аномального поведения кривой плавления Na.
Недавно мы опубликовали [6] результаты исследования плавления Na в ОЦК- и ГЦК-фазах, полученные на основе расчета динамики решетки методом линейного отклика ab initio в рамках теории функционала плотности (DFT) (описание метода см., например, в работе [7]). В рамках простого подхода, сочетающего DFT-расчеты (ab initio) фононных спектров с использованием эмпирического критерия Лиидемаиа, было дано количественное описание немонотонной кривой плавления Na. Было показано, что поведение Тт(р) хорошо объясняется изменением при сжатии фононного спектра твердой фазы без привлечения какой-либо информации о свойствах жидкой фазы. В частности, мы убедились, что наблюдаемое резкое понижение Тт(р) в
115
8*
ГЦК-структуро обусловлено смягчением при сжатии поперечной фононной моды и соответствующим возрастанием амплитуд тепловых колебаний атомов.
Наши DFT-расчеты Тт(р) натрия [6] были выполнены методом ah initio псевдопотеициала (AIPP), реализованном в программном комплексе Quantum Espresso [8], в квазигармоническом приближении, т. е. в предположении, что фононные частоты неявно зависят от температуры только за счет теплового расширения. В таком подходе не учитываются эффекты «истинного» энгармонизма, т.е. ангармонические взаимодействия фононных мод друг с другом и многофононные процессы. Эти эффекты могут приводить к существенной температурной зависимости фононных частот из-за конечного времени жизни фононов. Эффекты энгармонизма возрастают при увеличении температуры и пропорциональны среднеквадратичной амплитуде колебаний атомов. Высокие (около 1000 К) температуры исследуемой системы, как и смягченно поперечных фононов в ГЦК-фазе Na при сжатии, в принципе, могут приводить к заметным ангармоническим эффектам. Поэтому, несмотря на хорошее согласие наших квазигармонических расчетов [6] с экспериментом, оставался невыясненным вопрос о влиянии ангармонических взаимодействий на фазовую диаграмму Na в области высоких давлений и температур.
2. МЕТОДИКА РАСЧЕТА
Для исследования влияния энгармонизма на тер-модинэмические свойствэ натрия мы выполнили ah initio численное моделирование плавления Na при высоких давлениях с помощью квантовой молекулярной динамики, поскольку такой подход позволяет адекватно описывать свойствэ материалов при конечных температурах с полным учетом всех ангармонических эффектов. Результаты были частично представлены в краткой публикации [9]. Мы использовали метод молекулярной динамики Борнэ Оп-пенгеймерэ (BOMD) (см., например, [10]) в версии, реализованной в программном пэкете ('I'MI) [11]. В данном подходе движение атомов считается классическим и описывается обычными ньютоновскими уравнениями движения. Однако в отличие от стандартной молекулярной динамики с заданным априори феноменологическим потенциалом, в методе НОЛИ) взаимодействие между атомами рассчитывается ah initio непосредственно в процессе моделирования. Для этого путем минимизации функционала
электронной плотности вычисляются энергия электронной подсистемы при произвольном расположении попов, э также силы, действующие на ноны со стороны электронов. Данная процедура выполняется при каждом изменении положений попов. Такой подход, в отличие от широко распространенного метода Кара Парринелло, позволяет с хорошей точностью выполнять молекулярно-дннэмнческое моделирование ah initio не только для диэлектриков, но и в случае метэллов. Реализация данного подходэ требует очень больших вычислительных ресурсов, что накладывает серьезные ограничения на время моделирования и на возможное число частиц в исследуемой системе.
Интегрирование уравнений Ныотона в системе из N частиц с периодическими граничными условиями выполнялось с помогцыо разностной схемы Верле [12], которая дэет текущие координаты и скорости ионов с достаточно высокой точностью. В гэ-мильтоновой системе должна сохраняться полная энергия Е, э в состоянии равновесия также и температура Т. Использование конечно-рэзностного алгоритма интегрирования уравнений Ныотона приводит к возникновению флуктуэций Е, Т и других термодинамических величин вокруг их средних значений. Для более точного определения средних величин желательно использовать большое число частиц N в системе и малый временной шэг At. Таким образом, с учетом необходимости больших вычислительных зэтрэт, требуется тщательная оптимизация пэрэметров моделирования.
При моделировании ОЦК- и ГЦК-фэз Мэ были выбраны сверхъячейки, содержащие соответственно 4x4x4 (N = 128 этомов) и 3x3x3 (108 этомов) кубических элементарных ячеек. При расчете динамических корреляций, где необходима более высокая точность, для ГЦК-фазы использовалась сверхъячейка 4 x 4 x 4 (256 этомов). Временной шэг At выбирался равным 3.6 фс для О ЦК-фазы и 2.4 фс для ГЦК-фазы. Эти величины примерно отвечают значению /„,,„/•''>(). где /,„,„ минимальный период колебаний в исследуемой системе, который оценивается из максимальной частоты фононов. В практике подобных расчетов приемлемым обычно считается выбор _!///,„,„ ~ 0.01 0.1. Тестовые вычисления показали, что интересующие нас стэтические и динамические характеристики не изменяются при уменьшении шагэ до 1 фс. При каждом объеме моделирование проводилось для нескольких значений температуры, последовательно повышаемой вплоть до плавления системы. При этом, как правило, оказывалось достаточно выполнить около 6000 шэгов At.
что отвечает времени моделирования порядка 20 пс.
Минимизация функционала электронной плотности на каждом временном шаге проводилась методом ab initio псевдопотенциала в базисе плоских волн. При этом использовался сохраняющий норму нелокальный псевдопотенциал Труллера - Мартин-са с конфигурацией валентных электронов s1 и с учетом взаимодействия валентных и остовных электронов [11]. Обменно-корреляционный потенциал рассматривался в приближении локальной электронной плотности (LDA). Энергия обрезания плоских волн Ecut составляла 20 Ry. При повышении Ecut до 40 Ry траектории частиц, описывающие эволюцию системы, практически не изменялись. Интегрирование в обратном пространстве проводилось с использованием только Г-точки. Как показали наши тестовые расчеты, а также вычисления в работе [3], получаемые результаты в пределах численной точности совпадают с полученными при интегрировании по 2x2x2 = 8 особым k-точкам Монкхорста-Пака. Моделирование проводилось в iWT-ансамбле (N -число частиц, V — объем, Т — температура). Для термализации системы и изменения температуры в ней использовалась цепочка из четырех термостатов Нозе-Гувера [13,14] с временем релаксации около 0.1 пс.
20 t, пс
Рис.1. Среднеквадратичное отклонение атомов и(£) (пунктир) для ОЦК-натрия при р = 64.5 ГПа. Цифрами под графиком указаны температуры Т. Сплошной линией представлены значения и, полученные нами в квазигармоническом расчете [6] при тех же температурах. Стрелкой отмечено начало возрастания и{€) с увеличением времени моделирования, что отвечает возникновению диффузии в расплаве. Плавление системы происходит в интервале между 800 К и 900 К
3. РАСЧЕТ КРИВОМ ПЛАВЛЕНИЯ ^ МЕТОДОМ ВОМБ
Наиболее обоснованными подходами к определению температуры плавления Тш являются двухфазное моделирование или же сравнение вычисленных свободных энергий твердой и жидкой фаз. Однако эти методы слишком затратны, поэтому мы определяли наступление плавления более простым способом — по началу линейного роста среднеквадратичного смещения атомов (и2(Ь)) как функции времени моделирования Наши расчеты, как и имеющиеся в литературе данные, показывают, что так
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.