ЖЭТФ, 2019, том 155, вып. 5, стр. 878-885

СРАВНЕНИЕ РАЗЛИЧНЫХ МЕТОДОВ АТОМИСТИЧЕСКОГО

МОДЕЛИРОВАНИЯ ДЛЯ РАСЧЕТА ТЕМПЕРАТУРЫ

ФАЗОВОГО ПЕРЕХОДА НА ПРИМЕРЕ ЦИРКОНИЯ

И. С. Гордеевa,b*, С. В. Стариков^b,c

^a Московский физико-технический институт (государственный университет)

141701, Долгопрудный, Московская обл., Россия

^b Объединенный институт высоких температур Российской академии наук

125412, Москва, Россия

^c Interdisciplinary Centre for Advanced Materials Simulation, Ruhr-University

44801, Bochum, Germany

Поступила в редакцию 25 сентября 2018 г.,

после переработки 1 декабря 2018 г.

Принята к публикации 5 декабря 2018 г.

Выполнено сравнение различных методов расчета температуры фазового перехода на основе атоми-

стического моделирования, таких как двухфазное моделирование, однофазное моделирование и расчет

термодинамических потенциалов с помощью фононных спектров. На основе моделирования циркония

исследуются плавление и переход между различными кристаллическими фазами. Показано, что крите-

рий Линдемана является достаточно строгим при описании кривой плавления на фазовой диаграмме.

Однако пороговое значение параметра Линдемана, соответствующее фазовому переходу, может разли-

чаться в разных моделях. Кроме того, результаты исследования позволяют оценить точность расчета

температуры фазового перехода на основе термодинамического подхода с привлечением концепции эн-

тропии.

DOI: 10.1134/S0044451019050110

Метод молекулярной динамики (МД) [3] явля-

ется удобным средством, позволяющим исследовать

термодинамические эффекты в металлах на атомар-

1. ВВЕДЕНИЕ

ном уровне [4,5]. Важнейшим фактором в МД явля-

ется потенциал межатомного взаимодействия, кото-

Теоретические и экспериментальные исследова-

рый определяет силы взаимодействия [6-8]. Однако

ния фазовых диаграмм металлов являются одной

различные методы определения фазовых переходов

из основных целей материаловедения и атомисти-

могут давать различные результаты. Так, методом

ческого моделирования. Вопрос о построении фазо-

термодинамического интегрирования в работе [9]

вой диаграммы актуален, в частности, для цирко-

рассчитывается температура перехода между ОЦК-

ния. Данный металл широко используется в качест-

и ГПУ-фазами, равная 1233 K для потенциала EAM

ве конструкционного материала в атомной энерге-

[9], в то время как в работе [10] методом однофазно-

тике, медицине и аэрокосмической отрасли. Цирко-

го моделирования достигается значение 717 К для

ний имеет три твердотельные фазы: низкотемпера-

этого же потенциала. Кроме того, методом Монте

турную α-фазу с кристаллической ГПУ-решеткой;

Карло была получена [11] величина 1275 К. В рабо-

высокотемпературную β-фазу с кристаллической

те [12] метод МД используется для расчета порого-

ОЦК-решеткой; гексагональную ω-фазу, наблюдае-

вого значения параметра Линдемана, соответству-

мую при высоком давлении. Переход между ОЦК- и

ющего плавлению леннард-джонсовского кристалла

ГПУ-фазами циркония происходит при температуре

с открытой поверхностью. Рассчитанное пороговое

1130 K [1], плавление ОЦК-фазы — при 2125 K [2].

значение соответствует 0.22, что хорошо согласует-

ся с классической теорией Линдемана [13]. С дру-

^* E-mail: ilya.gordeyev@phystech.edu

878

ЖЭТФ, том 155, вып. 5, 2019

Сравнение различных методов атомистического моделирования.. .

гой стороны, в работе [14] квантовые расчеты фо-

2.1. Плавление

нонных свойств ОЦК-фазы циркония предсказыва-

ют довольно высокое значение параметра Линдема-

Двухфазное моделирование плавления являет-

на (около 0.35) при температуре плавления. Стоит

ся довольно строгим методом расчета температуры

отметить, что актуальным также является вопрос о

плавления. Основной идеей данного метода являет-

зависимости порогового значения данного парамет-

ся наблюдение за передвижением фазовой границы

ра от давления. Классическая теория предсказывает

двухфазной системы [20]. В расчетной ячейке созда-

неизменность этого параметра вдоль кривой плавле-

валась двухфазная система порядка 40000 атомов с

ния на фазовой диаграмме, однако сравнение с экс-

размерами 144 × 36 × 36Å3. После этого проводи-

периментальными данными [15, 16] показывает, что

лись наблюдения за движением фазовой границы.

это выполняется далеко не всегда.

Неподвижность ее свидетельствует о термодинами-

ческом равновесии, а движение — о преобладании

В настоящей работе рассчитываются температу-

одной фазы над другой. Длительность расчетов со-

ры переходов ОЦК-ГПУ и ОЦК-расплав методами

ставляла 150 пс. Варьируя температуры и давления

двухфазного и однофазного моделирования с помо-

и таким образом набирая статистику, мы получили

щью концепции энтропии на основании плотности

кривые плавления ОЦК-фазы (рис. 1). Температу-

фононных состояний. К наиболее обоснованным ме-

ры плавления при нулевом давлении для ADP-по-

тодам относят двухфазное и однофазное моделиро-

тенциала [17] составляет T_m = 1780 ± 30 K, а для

вание, которые непосредственно воспроизводят пе-

EAM-потенциала [9] — T_m = 2100 ± 30 K.

реход между фазами. Тем не менее данные мето-

Помимо этого была проведена проверка теории

ды чрезвычайно ресурсо- и времязатратные. Имен-

но эта отрицательная черта привела к созданию кос-

Линдемана. При классическом рассмотрении плав-

√

ление начинается, когда параметр X =

〈u²〉/W

венных методов расчета точек фазового перехода,

достигает значения X_m = 0.20-0.25 для большин-

таких как получение термодинамических потенциа-

ства твердых тел [14] (〈u²〉 — среднеквадратичное

лов (энтропии) через фононные спектры непосред-

смещение атома, W — средний радиус ячейки Виг-

ственно из МД. Кроме того, сам расчет температуры

нера - Зейтца для кристаллической решетки плавя-

перехода через термодинамические функции явля-

щейся фазы).

ется более привычным с точки зрения статистиче-

ской физики. Можно отметить, что прямые методы

моделирования переходов (двухфазное и однофаз-

T, K

ное МД-моделирование) являются вычислительны-

ми экспериментами, тогда как вычисление, основан-

ное на термодинамических функциях, больше соот-

2200

ветствует понятию «теоретический расчет». Таким

образом, результаты данной работы позволяют оце-

нить точность расчета энтропии на основе фонон-

ных свойств при описании фазового перехода.

2000

2. МЕТОДЫ И РЕЗУЛЬТАТЫ

В данной работе используется несколько методов

1800

атомистического моделирования для расчета тем-

ператур фазовых переходов в цирконии. Все рас-

четы были проведены с двумя различными потен-

циалами межатомного взаимодействия: EAM-потен-

P, ГПа

циалом из работы [9] (метод погруженного атома) и

Рис. 1. Кривые плавления циркония: 1 и 2 — результа-

ADP-потенциалом из работы [17] (потенциал с угло-

ты МД-расчета с потенциалами соответственно ADP [17]

вой зависимостью). Расчеты проводились с исполь-

и EAM [9]; 3 — температура плавления, измеренная экспе-

зованием кода LAMMPS [18], позволяющего иссле-

риментально [2]

довать широкий спектр задач [19].

879

И. С. Гордеев, С. В. Стариков

ЖЭТФ, том 155, вып. 5, 2019

По классической теории Линдемана [13] средне-

0.45

квадратичное смещение атомов задается в виде

9ℏ²T

〈u²〉 =

(1)

0.40

Mk_BΘ²

где M — масса одного атома, Θ — температура Де-

0.35

бая. Формулу (1) можно переписать для параметра

√

Линдемана

9ℏ²T

0.30

X =

(2)

W Mk_BΘ²

0.25

Подставляя в эту формулу дебаевскую темпера-

туру циркония Θ = 253 K, W = 1.77Å, получаем

теоретическую зависимость, которая при значении

0.20

X_m = 0.225 дает температуру 2125 К. Кроме того, из

EAM

теории Линдемана следует, что вдоль кривой плав-

0.15

ADP

T_m

ления на фазовой диаграмме параметр X_m остается

постоянным.

0.10

500

1000

1500

2000

2500

Рассчитывая с помощью кода LAMMPS сред-

T, K

неквадратичные смещения для обоих потенциалов,

мы получили зависимости параметра Линдемана от

0.50

температуры для ОЦК-фазы при различных давле-

ниях. Поскольку проведение двухфазного модели-

рования позволило нам рассчитать кривые плавле-

0.45

ния, мы смогли оценить значения параметра Линде-

мана при температуре плавления (рис. 2). Таким об-

разом, было получено соответствие прямого метода

0.40

вычисления кривой плавления (двухфазного моде-

лирования) с его косвенным аналогом (теория Лин-

демана). Причиной существенного завышения пара-

0.35

метра Линдемана в атомистическом моделировании,

вероятно, является тот факт, что плавление проис-

ходит из ОЦК-фазы, нестабильной при низкой тем-

0.30

пературе. Хотя температура Дебая для ГПУ-фазы

в комбинации с формулой (1) дает хорошее согласие

для температуры плавления, атомные колебания в

0.25

ОЦК-фазе значительно сильнее, чем допускает гар-

P, ГПа

моническое приближение.

Рис. 2. a) Зависимости параметра Линдемана X от темпе-

ратуры при P = 0: 1 и 2 — результаты расчетов для меж-

2.2. Фазовый переход между ОЦК- и

атомных потенциалов соответственно ADP и EAM; 3 —

ГПУ-фазами

теоретическая оценка, рассчитанная по формуле (2); се-

2.2.1. Двухфазное моделирование

рой областью отмечено характерное пороговое значение

при плавлении большинства металлов. б) Значения пара-

Для определения точки фазового перехода α-β

метра Линдемана X на рассчитанных кривых плавления в

используется вариация метода двухфазного модели-

МД-расчетах для потенциалов ADP (1) и EAM (2)

рования, которая предполагает, что переход осно-

ван на зародышеобразовании новой фазы. В расчет-

ной ячейке в наночастице одной фазы (70000 ато-

сматриваться в качестве катализатора для мартен-

мов, диаметр 15Å) создается сферическая область

ситного перехода (зародыш сферический, что суще-

другой фазы (3000 атомов, диаметр 4Å). Поскольку

ственно упрощает начало перехода). В процессе эво-

размеры зародыша новой фазы достаточно велики

люции двухфазная система переходит в однофаз-

по отношению к размерам системы (4 % от обще-

ную; длительность расчета не превышала 150 пс.

го объема), данная постановка задачи может рас-

Исчезновение зародыша свидетельствует о преоб-

880

ЖЭТФ, том 155, вып. 5, 2019

Сравнение различных методов атомистического моделирования.. .

Рис. 3. Рост зародыша ГПУ-фазы Zr для потенциала EAM [17] при температуре 590 K: черные атомы — ГПУ-фаза,

серые — ОЦК-фаза; a — t = 2.5 пс; б — t = 12.5 пс; в — t = 20 пс; г — t = 25 пс

ладании первоначальной фазы, а его увеличение,

ходится в метастабильных состояниях, что затруд-

в свою очередь, — об обратном (рис. 3). Визуали-

няет определение температуры фазового перехода.

зация проводилась с помощью программы OVITO

Однако, основываясь на том, что фазовый пере-

[21] (http://ovito.org/). Затем фазы меняются места-

ход начинается на границе зерен или на открытой

ми и проводятся аналогичные расчеты. В итоге по

поверхности, можно получить следующую форму-

температурным зависимостям (рис. 4) оценивают-

лу [10, 22]:

ся интервалы температур, из которых вычисляет-

∫

(

)

ся температура фазового перехода. Так, для потен-

(T -T₀)²

циала ADP [17] зародыш ГПУ-фазы начинает по-

exp

dT = CS,

(3)

k_BT

давляться при T > 1050 К, а зародыш ОЦК-фа-

зы — при T < 1350 К. Если принять за погреш-

где T_S — температура начала перехода при скорости

ность половину от разницы между указанными тем-

S нагрева/охлаждения, T₀ — равновесная темпера-

пературами, то температуру перехода можно оце-

тура фазового перехода, Q — энергия активации, ко-

нить как T^ADPα-β

= 1200 ± 80 K. Для потенциала

торая, как правило, принимает значение примерно

EAM [9] зародыш ГПУ-фазы начинает подавлять-

0.1 эВ, C — константа.

ся при T > 730 К, а зародыш ОЦК-фазы — при

Стоит отметить, что аналогичную формулу мож-

T < 1470 К. Следовательно, температуру перехода

но получить и для поверхностной диффузии. Кро-

можно оценить как T^EAMα-β = 1100 ± 150 K.

ме того, как грубую оценку эту формулу можно ис-

пользовать и для фазового перехода в бездефектном

2.2.2. Однофазное моделирование

кристалле. В последнем случае энергия активации

При охлаждении/нагреве в случае фазового пе-

является подгоночным параметром, характеризую-

рехода кристалл-кристалл однофазная система на- щим мобильность атомов на фазовой границе.

881

ЖЭТФ, вып. 5

И. С. Гордеев, С. В. Стариков

ЖЭТФ, том 155, вып. 5, 2019

T, K

E, эВ/атом

-6.26

1400

1300

-6.28

1200

-6.30

1100

-6.32

1000

350

400

450

500

550

100

t, пс

T, K

Рис. 4. Зависимости температуры системы от времени для

Рис. 5. Зависимости потенциальной энергии атомов от

потенциала ADP [17]: 1 — рост зародыша ОЦК-фазы; 2 —

температуры при охлаждении со скоростью S = 0.3 K/пс

подавление зародыша ОЦК-фазы; 3 — подавление заро-

при однофазном моделировании. Расчеты были выполне-

дыша ГПУ-фазы; 4 — рост зародыша ГПУ-фазы

ны с EAM-потенциалом [9]. Показаны различные зависи-

мости, различающиеся начальным микросостоянием мо-

делируемой системы

Основной рассчитываемой величиной в атоми-

Соответственно, температуру фазового перехода

стическом моделировании является зависимость

можно оценить как T^ADPα-β = 1195 ± 90 K и T^EAMα-β =

температуры перехода в охлаждаемой/нагреваемой

= 935 ± 240 K.

наночастице (с характерным размером около

10 нм) от скорости нагрева/охлаждения (рис. 5). В

2.2.3. Фононные спектры

качестве термостата для данных расчетов использо-

вался термостат Ланжевена. При проведении серий

Данный метод базируется на вычислении разно-

расчетов полученные температуры усредняются.

сти свободных энергий двух фаз путем расчета ста-

Соответственно, равновесная температура перехода

тистической суммы для фононов в гармоническом

рассчитывается как оптимальная температура,

приближении [23,24]. Свободную энергию F можно

обеспечивающая выполнение закона (3). Резуль-

записать как

таты аппроксимаций данных МД-моделирования

∫

[

(

)

представлены на рис. 6.

ℏω

F =

k_BT ln exp

2k_BT

Важно отметить, что в силу статистической при-

роды результатов МД-моделирования [5] температу-

(

)]

ℏω

ра перехода различается в расчетах даже с идентич-

- exp

g(ω) dω.

(4)

2k_BT

ными макропараметрами (начальная температура и

скорость нагрева/охлаждения). Поэтому необходи-

Если не рассматривать энергию нулевых колебаний,

мо проводить усреднение полученных температур

то энергия возбужденной системы фононов можно

перехода по начальным микросостояниям, что и бы-

записать в виде

ло проделано в данной работе. Для каждой скорос-

∫

ти S было проведено около 10 различных расчетов

ℏω

(

)

g(ω) dω,

(5)

и вычислялась средняя температура перехода.

ℏω

exp

-1

k_BT

Таким образом, были получены следующие ин-

тервалы температур: T^ADPα-β = 1020-1370 K для по-

причем плотность фононных состояний g(ω) пред-

тенциала ADP [17], T^EAMα-β = 460-1410 K для потен-

ставляет собой фурье-преобразование автокорреля-

циала EAM [9].

ционной функции [25]

882

ЖЭТФ, том 155, вып. 5, 2019

Сравнение различных методов атомистического моделирования.. .

T_S, K

500

500 K

1550

460

450

1450 K

1500

1410

400

430

1450

1390

1400

350

0.5

1.0

1.5

0.5

1.0

1.5

S, K/пс

T_S, K

1050

1060 K

1020

1500

1400 K

1000

1370

990

1330

1400

950

900

1300

0.5

1.0

1.5

0.5

1.0

1.5

S, K/пс

Рис. 6. Зависимости наблюдаемой температуры перехода TS от скорости охлаждения/нагрева S: символы соответствуют

результатам однофазного моделирования; линии — аппроксимациям по формуле (3) с различными значениями парамет-

ра T0 (приведены на рисунках). Показано несколько типов расчетов: а — охлаждение ОЦК-фазы для потенциала EAM

[9]; б — нагрев ГПУ-фазы для потенциала EAM [9]; в — охлаждение ОЦК-фазы для потенциала ADP [17]; г — нагрев

ГПУ-фазы для потенциала ADP [17]

∫T

v_j — скорость j-го атома. Плотность фононных со-

g(ω) = 2

VACF(t) cos(ωt) dt,

(6)

стояний рассчитывается с помощью кода LAMMPS;

усреднение проводится по всем атомам в ячейке.

Вид рассчитанного спектра для ОЦК- и ГПУ-фаз

где

представлен на рис. 7.

Важно отметить, что расчет термодинамических

∑

функций на основе фононных спектров (как и рас-

〈v_j (t)|v_j (0)〉

чет самих этих спектров) базируется на гармони-

VACF(t) =j=1

(7)

ческом приближении в описании колебаний атомов.

∑

В то же время, цирконий обладает сильной ангар-

〈v_j (0)|v_j (0)〉

моничностью атомных колебаний, что существен-

j=1

883

И. С. Гордеев, С. В. Стариков

ЖЭТФ, том 155, вып. 5, 2019

g, отн. ед.

Таблица. Результаты моделирования фазового пе-

рехода между фазами α и β в цирконии

0.25

Потенциал

0.20

Метод расчета

EAM

ADP

T_α-β, K ΔT, K T_α-β, K ΔT, K

0.15

Двухфазный

1100

150

1200

Однофазный

935

240

1195

0.10

Фононы

1180

1520

0.05

ΔF = ΔE_MD - T ΔS_phonon,

(8)

где

, ТГц

Рис. 7. Плотности фононных состояний g(ω) при темпе-

ΔS_phonon =

(ΔF_phonon - ΔE_phonon) ,

(9)

ратуре T = 700 K в случае потенциала EAM [9] для фаз

ОЦК (1) и ГПУ (2)

ΔE_MD — разность потенциальных энергий фаз,

F, эВ

рассчитанная непосредственно в ходе МД-модели-

0.02

рования, а ΔE_phonon, ΔF_phonon и ΔS_phonon — соот-

ветственно разности энергий, свободных энергий и

энтропий, рассчитанные через плотность фононных

0.01

состояний.

В работе [28] на примере фазовых переходов в ти-

тане было показано, что модели, основанные только

на гармонических колебаниях, дают сильно завы-

шенную температуру перехода ГПУ-ОЦК. В насто-

ящей работе, как и в работе [28], энергия системы

-0.01

берется напрямую из МД-расчета, что существенно

увеличивает точность расчета температуры перехо-

-0.02

да. Отметим, что при описании взаимодействия ато-

мов с помощью потенциала ADP [17], они обладают

большим ангармонизмом колебаний в сравнении со

-0.03

случаем потенциала EAM [9].

Расчеты проводились для 32000 атомов в тече-

500

1000

1500

ние 240 пс. В результате были получены зависи-

T, K

мости разности свободных энергий от температуры

Рис. 8. Рассчитанные разности свободных энергий фаз

(рис. 8) и сделана оценка температуры фазового пе-

ОЦК и ГПУ, ΔF = FHCP - FBCC , для потенциалов EAM

рехода между фазами α и β для обоих используе-

(1) и ADP (2)

мых потенциалов: T^ADPα-β = 1520 ± 20 K, T^EAMα-β =

= 1180 ± 20 K.

но при высоких температурах [26]. Следствием это-

го является существенная зависимость теплоемко-

3. ОБСУЖДЕНИЕ

сти от температуры [27]. Учет ангармонических эф-

фектов в формулах (4) и (5) является нетривиаль-

В настоящей работе получены оценки температу-

ной задачей, вследствие чего появляется необходи-

ры фазового перехода между ОЦК- и ГПУ-фазами

мость ввода поправки на ангармонизм. В данной ра-

в цирконии для двух потенциалов, EAM [9] и ADP

боте разность между свободными энергиями ОЦК-

[17] (таблица). Результаты демонстрируют, что рас-

и ГПУ-фаз рассчитывалась по формуле:

чет термодинамических потенциалов, основанный

884

ЖЭТФ, том 155, вып. 5, 2019

Сравнение различных методов атомистического моделирования.. .

на вычислении фононных спектров и концепции эн-

J. J. Möller, M. Mrovec, I. Bleskov et al., Phys. Rev.

тропии, при слабом ангармонизме потенциала взаи-

Mater. 2, 093606 (2018).

модействия (потенциал EAM [9]) хорошо согласует-

A. Rohskopf, H. R. Seyf, K. Gordiz et al., NPJ

ся с прямыми методами, такими как двухфазное и

Comput. Mater. 3:27 (2017).

однофазное моделирования. В обратном случае (по-

тенциал ADP [17]) температура перехода получает-

M. I. Mendelev and G. J. Ackland, Phil. Mag. Lett.

87, 349 (2007).

ся завышенной.

Стоит отметить и тот факт, что методы фонон-

10.

R. W. Smith, Comput. Condens. Matter 13,

ных спектров и двухфазного моделирования в дан-

(2017).

ной интерпретации значительно менее ресурсоем-

кие, чем метод однофазного моделирования. Одна-

11.

T. L. Underwood and G. J. Ackland, Comp. Phys.

ко, основываясь на плотности фононных состояний,

Comm. 215, 204 (2017).

можно проводить расчеты точек фазового перехода

12.

Z. H. Jin, P. Gumbsch, K. Lu, and E. Ma, Phys. Rev.

при различных давлениях. С другой стороны, двух-

Lett. 87, 055703 (2001).

фазное/однофазное моделирование фазового пере-

хода между кристаллическими фазами требует на-

13.

F. A. Lindemann, Physik. Z. 11, 609 (1910).

личия открытой поверхности для релаксации объе-

14.

I. Schnell and R. C. Albers, J. Phys.: Condens. Matter

ма (или формы) и поэтому осуществимо только для

18, 1483 (2006).

нулевого давления.

Также в работе с помощью метода двухфазного

15.

G. H. Wolf and R. Jeanloz, J. Geophys. Res. 89, 7821

моделирования были рассчитаны кривые плавления

(1984).

циркония для двух разных потенциалов. Показано,

16.

L. Koci, R. Ahuja, L. Vitos et al., Phys. Rev. B 77,

что параметр Линдемана не меняется вдоль кривой

132101 (2008).

плавления, однако само пороговое значение при тем-

пературе плавления сильно различается для разных

17.

D. E. Smirnova and S. V. Starikov, Comput. Mate-

rials Sci. 129, 259 (2017).

потенциалов.

Расчеты проводились с использованием вычис-

18.

S. J. Plimpton, J. Comp. Phys. 117, 1 (1995).

лительного кластера МВС-10П (Межведомствен-

19.

А. В. Янилкин, П. А. Жиляев, А. Ю. Куксин и др.,

ный суперкомпьютерный центр РАН).

Вычислит. методы и программир. 11, 111 (2010).

Финансирование. Работа выполнена при под-

20.

S. V. Starikov and V. V. Stegailov, Phys. Rev. B 80,

держке по Программе финансирования исследова-

220104(R) (2009).

ний Президиума Российской академии наук № I.13

21.

A. Stukowski and K. Albe, Modelling Simul. Mater.

(координатор академик В. Е. Фортов).

Sci. Eng. 18, 085001 (2010).

22.

Y. T. Zhu and J. H. Devletian, Metallurgical Trans-

ЛИТЕРАТУРА

actions A 22, 1993 (1991).

1. V. A. Petukhov, High Temp.-High Press. 35/36, 15

23.

W. Cochran and R. A. Cowley, Phonons in Per-

(2003).

fect Crystals, S. Flogge Encyclopedia of Physics,

Vol. XXV/2a (1967), pp. 59-156.

2. D. A. Young, Phase Diagrams of the Elements, Law-

rence Livermore Laboratory (1975).

24.

Л. Д. Ландау, Е. М. Лифшиц, Статистическая

физика. Часть 1, Физматлит, Москва (2010).

3. Г. Э. Норман, В. В. Стегайлов, Матем. моделир.

24, 3 (2012).

25.

D. V. Minakov, P. R. Levashov, and V. B. Fokin,

Comput. Materials Sci. 127, 42 (2017).

4. С. В. Дмитриев, Е. А. Корзникова, А. П. Четвери-

26.

T. May, W. Müller , and D. Strauch, Phys. Rev. B 57,

ков, ЖЭТФ 153, 417 (2018).

5758 (1998).

5. Г. Э. Норман, С. В. Стариков, В. В. Стегайлов,

27.

D. E. Smirnova, S. V. Starikov, and I. S. Gordeev,

ЖЭТФ 141, 910 (2012).

Comput. Materials Sci. 152, 51 (2018).

6. S. Rawat and N. Mitra, Comput. Mater. Sci. 126,

28.

K. Masuda-Jindo, S. R. Nishitani, and Vu Van Hung,

228 (2017).

Phys. Rev. B 70, 184122 (2004).

885