О применимости потенциалов модели погруженного атома (ЕАМ) к жидким кремнию и германию
- Авторы: Белащенко Д.К.1
-
Учреждения:
- Национальный исследовательский технологический университет МИСИС
- Выпуск: Том 99, № 1 (2025)
- Страницы: 122-134
- Раздел: ХЕМОИНФОРМАТИКА И КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ
- Статья получена: 17.04.2025
- Статья одобрена: 17.04.2025
- Статья опубликована: 17.04.2025
- URL: https://ogarev-online.ru/0044-4537/article/view/288152
- DOI: https://doi.org/10.31857/S0044453725010122
- EDN: https://elibrary.ru/EIAKHV
- ID: 288152
Цитировать
Аннотация
Предложены потенциалы модели погруженного атома (ЕАМ – embedded atom model) для жидких кремния и германия. Потенциалы рассчитаны по дифракционным данным с помощью алгоритма Шоммерса и представлены в виде таблиц, а также в виде кусочно-непрерывных полиномов. Каждый парный вклад в потенциал имеет вид, близкий к жесткосферному со ступенькой вниз. Рассчитаны свойства жидких Si и Ge при температурах до 2000 К: плотность, энергия, модуль всестороннего сжатия, коэффициенты самодиффузии. Отмечено, что согласие с опытом хорошее. Установлено, что при обычных плотностях жидких Si и Ge направленность связи практически полностью исчезает после плавления. Предполагается, что направленность связи может появиться при нагревании и уменьшении плотности расплавов в 2–3 раза.
Ключевые слова
Полный текст
ВВЕДЕНИЕ
Ниже рассмотрена проблема компьютерного моделирования жидкостей, получающихся после плавления и частичной металлизации веществ с ковалентной связью. Такие простые вещества расположены в полосе Периодической системы, включающей кремний и германий (4-я группа), сурьму (5-я группа) и теллур (6-я группа). Большинство этих веществ, в том числе кремний и германий являются в твердом состоянии полупроводниками, а в жидком состоянии проявляют металлические свойства. В случае сурьмы вклад ковалентной связи в кристалле невелик, а в теллуре ковалентность довольно значительна. При моделировании кристаллических ковалентных веществ обычно применяются потенциалы, стимулирующие направленность (наличие валентных углов). Однако в жидком состоянии признаки ориентационного взаимодействия наблюдаются не всегда. Поэтому возникает вопрос, можно ли применять один и тот же межчастичный потенциал при моделировании такой системы и в кристаллическом, и в жидком состоянии. Если типы связи в кристалле и в жидкости реально различны, то и потенциалы межчастичного взаимодействия должны отличаться. На это указывают, в частности, расчеты ab initio, из которых видны кардинальные различия структуры кристалла и жидкости.
В случае кремния и германия попытки применения единых потенциалов для кристалла и жидкости были довольно часты, однако хороших результатов на этом пути добиться не удалось. Ниже будут рассмотрены возможности моделировать указанные вещества в жидком состоянии с применением обычных сферически симметричных потенциалов, используемых для металлов.
МЕЖЧАСТИЧНЫЕ ПОТЕНЦИАЛЫ ЕАМ ДЛЯ КРЕМНИЯ И ГЕРМАНИЯ
Моделирование кремния. Твердый кремний – полупроводник. Кристаллическая решетка кремния кубическая с координационным числом КЧ = 4. Каждый атом находится в центре тетраэдра и связан с соседними атомами ковалентными связями типа sp3. Длина этой связи равна 2.352 Å. Дифракционные данные о структуре жидкого кремния приведены в книге Й. Васеды [1], а при давлениях 4–23 ГПа – в [2, 3]. Парная корреляционная функция (ПКФ) кремния для температуры 1733 К показана на рис. 1. Первый пик довольно высок и узок, так что первая координационная сфера хорошо определена; ее радиус равен 3.05 Å. Среднее координационное число КЧ = 6.24 ± 1.45. Имеются данные о плотности твердого и жидкого кремния до 2000 К [4, 5], об энергии кремния до 3600 К [6, 7], сжимаемости [8, 9], коэффициентах самодиффузии [3, 10], поведении при высоких давлениях [3, 11, 12] и др.
Рис. 1. Парная корреляционная функция Si при 1773 К: 1 – дифракционные данные [1], 2 –модель с потенциалом ЕАМ.
Для моделирования веществ с решеткой алмаза – углерода, кремния, германия – методом молекулярной динамики (МД) был предложен ряд межчастичных потенциалов, которые реализуют направленность связи: Стиллинжера–Вебера (SW) [13–15], Терсофа (Ter) [16–18] и другие, зависящие от локальной структуры, например “environment dependent interatomic potential” (EDIP) [12, 19, 20]. Такие потенциалы включают валентные углы между ближайшими 4 соседями центрального атома, что позволяет реализовать тетраэдрическую конфигурацию. Значительное число потенциалов указано в работе [21] и приведено в репозитории NIST [22]. Потенциалы, учитывающие валентные углы, обычно используют для расчетов структуры, коэффициентов диффузии и температуры плавления, но не энергетических характеристик. М. Баскес и др. [23, 24] применили для неизотропных веществ потенциал МЕАМ (modified embedded atom model). Довольно широко применяется метод ab initio [17, 21]. Данные, полученные этим методом, использовали для расчета потенциалов программой типа Potfit [25, 26]. В упомянутых работах были исследованы преимущественно модели кристаллического кремния и его сплавов.
Переход к анализу жидкого кремния вызывает трудности из-за существенного изменения типа связи. В случае потенциалов SW и Ter максимальное число атомов модели имеет КЧ = 6 (44% атомов в случае SW и 33% в случае Ter), а число атомов с КЧ = 4 невелико – 7% в SW и 8% в Ter, см. рис. 3 в [17]). Расчетные температуры плавления кремния с этими потенциалами составляют, соответственно, 1400 и 2750 К (фактическое значение 1687 К). В этих условиях применение потенциалов типа SW и Ter сомнительно. Что касается потенциала МЕАМ, то он для жидкого кремния в [23, 24] не применялся.
Рис. 2. Парный вклад в потенциал ЕАМ жидкого кремния. Алгоритм Шоммерса.
Рис. 3. Частота появления КЧ Z в модели жидкого кремния при 1690 К. Радиус координационной сферы 3.05 Å: 1 – потенциал ЕАМ, 2–1800 K, метод ab initio [21].
В связи с этом в настоящей работе предложен потенциал модели погруженного атома – ЕАМ (embedded atom model [27]) для жидкого кремния с использованием дифракционных данных. Потенциальная энергия металла записывается в виде:
(1)
Здесь Φ(ρi) – потенциал погружения i-го атома, зависящий от “эффективной электронной плотности” ρ в месте нахождения центра атома, а вторая сумма по парам атомов содержит обычный парный потенциал. Эффективная электронная плотность в точке нахождения атома создается окружающими атомами и определяется по формуле:
(2)
где ψ(rij) – вклад в электронную плотность атома i от соседа номер j. В расчетах используются три подгоночные функции Φ(ρ), ϕ(r) и ψ(r), так что возможности согласования расчетных свойств с экспериментальными значениями очень широки.
Для нахождения парного вклада ϕ(r) применим алгоритм Шоммерса [28, 29], в котором поправка к текущей версии парного вклада в потенциал определяется разницей между дифракционной и модельной ПКФ. Ниже будем определять расстояние между двумя графиками – гистограммами ПКФ g0(ri) и g(ri) – через невязку Rg:
(3)
где g0 (rj) – гистограмма целевой ПКФ, g (r) – гистограмма ПКФ модели, n1 и n2 – границы суммирования табличных данных, а j – номер элемента гистограммы ПКФ. Вначале с помощью ПКФ Й. Васеды [1] была построена алгоритмом Шоммерса модель жидкого кремния при 1733 К c парным потенциалом и очень небольшой невязкой Rg = 0.023 (см. рис. 1). Полученный при этом парный вклад в потенциал ЕАМ приведен в табл. 1 и на рис. 2 с шагом dr = 0.05 Å. Он похож на жесткосферный со ступенькой вниз. Радиус обрыва rc = 8.855 Å.
Таблица 1. Парный вклад в потенциал жидкого кремния. Алгоритм Шоммерса при 1733 К. Парная корреляционная функция [1], невязка Rg = 0.023
r, Å | ϕ(r), эВ | r, Å | ϕ(r), эВ | r, Å | ϕ(r), эВ | r, Å | ϕ(r), эВ |
1.5 | 122.291000 | 3.35 | –0.073576 | 5.2 | –0.000563 | 7.05 | –0.000178 |
1.55 | 101.294000 | 3.4 | –0.073472 | 5.25 | –0.002079 | 7.1 | –0.000617 |
1.6 | 82.287700 | 3.45 | –0.072458 | 5.3 | –0.003564 | 7.15 | –0.000834 |
1.65 | 65.261100 | 3.5 | –0.071103 | 5.35 | –0.005022 | 7.2 | –0.000750 |
1.7 | 50.214400 | 3.55 | –0.069557 | 5.4 | –0.006002 | 7.25 | –0.000432 |
1.75 | 37.180300 | 3.6 | –0.067373 | 5.45 | –0.006646 | 7.3 | 0.000022 |
1.8 | 26.239600 | 3.65 | –0.064463 | 5.5 | –0.007191 | 7.35 | 0.000467 |
1.85 | 17.473600 | 3.7 | –0.060170 | 5.55 | –0.007183 | 7.4 | 0.001153 |
1.9 | 10.839500 | 3.75 | –0.054712 | 5.6 | –0.006751 | 7.45 | 0.001830 |
1.95 | 6.067210 | 3.8 | –0.048719 | 5.65 | –0.005952 | 7.5 | 0.002192 |
2.00 | 2.767880 | 3.85 | –0.042092 | 5.7 | –0.005362 | 7.55 | 0.002491 |
2.05 | 0.755346 | 3.9 | –0.036168 | 5.75 | –0.005026 | 7.6 | 0.002586 |
2.1 | 0.159373 | 3.95 | –0.030396 | 5.8 | –0.004385 | 7.65 | 0.002628 |
2.15 | –0.001554 | 4.00 | –0.023942 | 5.85 | –0.003832 | 7.7 | 0.002827 |
2.2 | –0.101008 | 4.05 | –0.016892 | 5.9 | –0.003009 | 7.75 | 0.002638 |
2.25 | –0.152195 | 4.1 | –0.009745 | 5.95 | –0.002833 | 7.8 | 0.002702 |
2.3 | –0.189574 | 4.15 | –0.002723 | 6 | –0.002206 | 7.85 | 0.002747 |
2.35 | –0.213938 | 4.2 | 0.003849 | 6.05 | –0.001500 | 7.9 | 0.002539 |
2.4 | –0.228382 | 4.25 | 0.009620 | 6.1 | –0.000687 | 7.95 | 0.002463 |
2.45 | –0.234435 | 4.3 | 0.014813 | 6.15 | –0.000636 | 8 | 0.002386 |
2.5 | –0.232160 | 4.35 | 0.019150 | 6.2 | –0.000856 | 8.05 | 0.002338 |
2.55 | –0.224053 | 4.4 | 0.023231 | 6.25 | –0.001311 | 8.1 | 0.002013 |
2.6 | –0.209804 | 4.45 | 0.026716 | 6.3 | –0.001707 | 8.15 | 0.001976 |
2.65 | –0.191907 | 4.5 | 0.029711 | 6.35 | –0.002199 | 8.2 | 0.001765 |
2.7 | –0.171074 | 4.55 | 0.032243 | 6.4 | –0.002325 | 8.25 | 0.001668 |
2.75 | –0.150465 | 4.6 | 0.033630 | 6.45 | –0.002367 | 8.3 | 0.001353 |
2.8 | –0.129192 | 4.65 | 0.034171 | 6.5 | –0.001990 | 8.35 | 0.001162 |
2.85 | –0.108827 | 4.7 | 0.033699 | 6.55 | –0.001173 | 8.4 | 0.000972 |
2.9 | –0.094394 | 4.75 | 0.032001 | 6.6 | –0.000397 | 8.45 | 0.000613 |
2.95 | –0.086053 | 4.8 | 0.028947 | 6.65 | 0.000469 | 8.5 | 0.000282 |
3 | –0.081530 | 4.85 | 0.025344 | 6.7 | 0.001206 | 8.55 | 0.000034 |
3.05 | –0.077580 | 4.9 | 0.020481 | 6.75 | 0.001713 | 8.6 | –0.000053 |
3.1 | –0.074926 | 4.95 | 0.015896 | 6.8 | 0.002049 | 8.65 | –0.000177 |
3.15 | –0.073313 | 5 | 0.011360 | 6.85 | 0.001965 | 8.7 | –0.000182 |
3.2 | –0.073026 | 5.05 | 0.007204 | 6.9 | 0.001731 | 8.75 | –0.000497 |
3.25 | –0.072734 | 5.1 | 0.003855 | 6.95 | 0.000888 | 8.8 | –0.000248 |
3.3 | –0.073358 | 5.15 | 0.001455 | 7 | 0.000334 | 8.85 | 0.000000 |
Потенциал погружения имеет вид:
(4)
Коэффициенты a2 и b2 связаны с a1 и c1 условием непрерывности функции Φ(ρ) и ее производной: a2 = a1 + c1 (ρ1 – ρ0)2 и b2 = 2c1(ρ1 – ρ0). Величины a1, c1, c2, ρ1 должны быть заданы. Параметры потенциала погружения ЕАМ Φ(ρ) приведены в табл. 2. Значения ρ > ρ2 могут быть получены только в сильно сжатых состояниях.
Таблица 2. Жидкий кремний. Параметры потенциала погружения F(r)
p1 | 0.46785 | n | 2048 |
p2 | 0.80 | dr, Å | 0.05 |
r1 | 0.875 | r0, Å | 9.86 |
r2 | 1.125 | rc, Å | 8.855 |
a1, эВ | –3.2823 | c1, Å | 2.05 |
С помощью потенциала ЕАМ была построена серия моделей жидкого кремния размером 2048 атомов в основном кубе с периодическими граничными условиями, при температурах 1690–2000 К, в режимах NVT и NpT. Парная корреляционная функция g(r) при 1733 К показана на рис. 1 в отличном согласии с дифракционными данными [1].
Приведенное на рис. 3 распределение КЧ модели кремния имеет вид, характерный для плотных упаковок, и КЧ = 4 ничем не выделяется. Это распределение близко к полученному методом ab initio при 1800K [21]. Цепочечная или тетраэдрическая упаковка в жидком кремнии не видна. Плавление разрушает сетевую структуру.
Распределение азимутальных углов (рис. 4) также не похоже на случай потенциала Ter [17, 30], где виден острый максимум при углах ~60°. Расчеты методом ab initio [21] приводят к кривой с двумя максимумами при значениях углов 60° и 90°, как и в случае ЕАМ (рис. 4), однако правый максимум при 90° в [21] немного выше левого пика. Потенциалы SW и Ter приводят к угловым распределениям, сильно отличающимся по форме от кривой на рис. 4.
Рис. 4. Азимутальные углы θ в модели Si при 1733 К.
Результаты расчетов свойств жидкого кремния при температурах до 2000 К приведены в табл. 3. Как видно из таблицы, потенциал ЕАМ позволяет получить правильную зависимость плотности жидкого кремния от температуры до 2000 К, хорошее согласие ПКФ с опытом при 1733 К (рис. 1), хорошее согласие с опытом модуля всестороннего сжатия. При 1733 К теплоемкости Cv и Cp модели кремния равнялись 16.9 и 21.4 Дж/(моль К). Пересчет адиабатического модуля Ks = 35.7 ГПа [8, 9] на изотермический KT был проведен с учетом отношения теплоемкостей Cp/Cv = 1.266; в итоге КТ = 28.2 ГПа. При значении c1 = 2.05 эВ величина модуля, найденная по зависимости давления от объема, равна КТ = 30.6 ГПа. При нагревании от 1690 до 2000 К модуль КТ уменьшается довольно медленно – всего на 10%. Хорошее согласие с опытом [7] получено для энергии жидкого кремния. С ростом температуры наблюдается небольшое занижение энергии моделей, обусловленное неучетом электронной теплоемкости моделей кремния. Электронная теплоемкость при 1600–2000 К составляет ~2 Дж/(моль К).
Таблица 3. Свойства моделей Si, полученные методом МД при p ∼ 0.001 ГПа, p1= 0.46785, p2 = 0.800, c1= 2.0500
T, K | d, г/см3 | <ρ>b | Rg | –E, кДж/моль | KT, ГПа | D×105, см2/с | |||||
МД | Эксп. [4] | Эксп. [5] | –EMD | –Eexp [7] | МД | MD | |||||
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
0 | – | – | – | – | – | – | 445.67 [6] | – | – | – | – |
298a | 2.329 | – | 2.328 | 1.0153 | – | 358.30 | 442.46 | – | 98.74 | – | – |
1000a | – | – | 2.320 | – | – | – | 425.37 | – | – | – | – |
1300a | – | – | 2.305 | – | – | – | 417.22 | – | – | – | – |
1500a | – | – | 2.297 | – | – | – | 411.61 | – | – | – | – |
1690a | – | – | 2.285 | – | – | – | 405.25 | – | – | – | – |
1690 liq | 2.525 | – | – | 1.0185 | – | 355.90 | 356.00 | 30.4 | 28.2 | 25.3 | – |
1700 | 2.523 | 2.547 | – | 1.0211 | – | 355.71 | 355.65 | – | 37.4 | 24.2 | – |
1733 | 2.517 | 2.515 | 1.0145 | 0.030 | 354.92 | 354.75 | 30.6 | – | 26.5 | – | |
1750 | 2.516 | 2.534 | – | 1.0159 | – | 354.60 | 354.29 | – | 38.6 | – | – |
1800 | 2.508 | 2.520 | – | 1.0158 | – | 353.52 | 352.93 | 27.3 | 39.4 | 26.1 | 19–20 |
1900 | 2.491 | 2.494 | – | 1.0038 | – | 351.47 | 350.21 | 39.8 | 26.6 | – | |
2000 | 2.470 | 2.467 | – | 0.9999 | – | 349.34 | 347.49 | 27.6 | – | 32.3 | – |
a Кристалл.
b Стандартные отклонения растут сверху вниз от 0.073 до 0.079.
Коэффициент самодиффузии D находили по зависимости среднего квадрата смещения частиц от времени. Потенциал ЕАМ дает высокие значения коэффициента D ~2×10–4 см2/с (табл. 3). Метод ab initio дает значение D = 2.26×10–4 такого же порядка [21, 32]. Однако, потенциалы SW и Ter приводят к коэффициентам самодиффузии жидкого кремния порядка 10–5 см2/с [30, 33]. Эффект повышенного коэффициента самодиффузии (до ~10–4 см2/с) наблюдался ранее у моделей жидкой сурьмы при Т > 1000 K [34], но он при нагревании постепенно исчезает.
Близость потенциала кремния к жесткосферному со ступенькой вниз (рис. 2) позволяет проверить формулу для коэффициента самодиффузии модели жестких сфер [35]:
(5)
где n = (N/V) σ3, N – число частиц в объеме V, σ – диаметр сферы, m – ее масса, D0 = (3/8) σ(kBT/πm)1/2.Умножая числитель и знаменатель дроби в скобке на число Авогадро, получаем при T = 1733 К следующее: σ = 2.20 Å, N/V = 0.053932 ат/Å3, n = 0.5743, D0 = (3/8) 2.2×10–8×(8.31×107×1733/π/28.085)1/2 = 3.33×10–4 см2/с, и, наконец, DHS = 2.86×10–4 см2/с. Эта величина хорошо согласуется с цифрой 2.65×10–4, полученной методом МД, и со значением D, найденным методом ab initio (2.02·10–4 см2/с при 1800°С [21]).
Механизм самодиффузии можно уточнить, проверив выполнимость для жидкого кремния уравнения Стокса–Эйнштейна, связывающего самодиффузию с вязкостью:
, (6)
где η – динамическая вязкость и ra – радиус атома. Выбирая для вязкости при 1700 К величину 6.05 мПуаз [4] и ra = 1.1×10–8 см (см. выше), получаем D = 1.9×10–4 см2/с, что согласуется с фактическим значением (1.9–2.0)×10–4 (табл. 3). Выполнимость уравнения Стокса–Эйнштейна означает отсутствие признаков ассоциации реальной жидкости и, соответственно, компьютерной модели.
Итак, результаты моделирования показывают, что свойства жидкого кремния хорошо описываются в рамках теории простых жидкостей с потенциалом, близким к жесткосферному, с диаметром сфер 2.20 Å. Этот диаметр меньше межатомного расстояния в кристалле (2.352 Å).
При моделировании жидкого кремния можно использовать парный вклад в потенциал ЕАМ, заданный не в виде таблицы данных, а в виде кусочно-непрерывного полинома [29, 36]. Выделим на оси абсцисс r последовательные k интервалов точками деления r1 > 0, r2, r3 …rk+1. Координата r1 должна быть немного меньше минимального межатомного расстояния rmin. Первый интервал заключен между точками r1 и r2. Координата rk+1 должна быть равна радиусу обрыва потенциала rc. Алгоритм Шоммерса дает значения потенциала в интервалах от 1-го до k-го, то есть при r1 ≤ r ≤ rk+1, а на расстояниях r < r1 не работает. При r > r1 потенциал описывается выражением:
(7)
где i – номер интервала на оси расстояний, dim – коэффициенты разложения в ряд, и функция Хевисайда равна 1 в интервале ri ≤ r ≤ ri+1 и нулю в остальных интервалах. При всех r = ri непрерывен сам потенциал и его первая производная. Расстояние r задается в Å. Степень полиномов равна n. В случае кремния выбрано n = 7, k = 4. Коэффициенты dim получены в программе аппроксимации данных табл. 1 и приведены в табл. 4. В этой процедуре значения потенциала были надежно определены для расстояний r ≥ r1 = 2.15 Å, которые реально встречались в модели кремния при 1733 K.
Таблица 4. Коэффициенты разложения парного вклада в потенциал ЕАМ кремния
aim | Номер интервала i / Границы интервала ri – ri+1, Å | |||
1 / 2.10–2.90 | 2 / 2.90–4.20 | 3 / 4.20–6.00 | 4 / 6.00–8.85 | |
ai0 | –0.97113050520420D-01 | 0.31690669711679D-02 | –0.31873278785497D-02 | 0.00000000000000D+00 |
ai1 | 0.21899378299713D+00 | 0.12773799896240D+00 | 0.34430783707649D-02 | 0.00000000000000D+00 |
ai2 | –0.62291081298754D+01 | 0.27154179420195D-01 | 0.84230139063026D-01 | –0.35165857912545D-01 |
ai3 | –0.79015407873123D+02 | 0.43481827129286D+00 | 0.56141844748633D+00 | –0.13246413781381D+00 |
ai4 | –0.41364497186011D+03 | 0.13526037867501D+01 | 0.11214214299747D+01 | –0.16534787554672D+00 |
ai5 | –0.10115477873406D+04 | 0.14864527900288D+01 | 0.89595677539435D+00 | –0.93427224218695D-01 |
ai6 | –0.11462934613386D+04 | 0.74986992706763D+00 | 0.31125775931700D+00 | –0.24739723425154D-01 |
ai7 | –0.48846592560798D+03 | 0.15751335772576D+00 | 0.39450678719440D-01 | –0.24982008692301D-02 |
Поскольку при высоких температурах и давлениях минимальные межатомные расстояния rmin убывают, то необходимо провести экстраполяцию потенциала на расстояния, меньшие rmin. Для кремния принято, что при r < r1:
где r выражено в Å и ϕ(2.15 Å) = –0.0787061 эВ. Переход к записи потенциала в форме (7) приводит к небольшой потере точности, так что невязка Rg функции g(r) кремния при 1733 К увеличивается от 0.023 до ~0.06.
Моделирование германия. Твердый германий – полупроводник. Кристаллическая решетка такая же, как у кремния. Дифракционные данные о структуре жидкого германия приведены в работах [1, 33, 37], а при давлениях до 25 ГПа – в [38]. Парная корреляционная функция германия при температуре 1253 К показана на рис. 5. Первый пик при 2.75 Å довольно высок и узок, так что первая координационная сфера хорошо определена с радиусом 3.60 Å. Среднее координационное число КЧ = 8.63 ± 1.41. Имеются данные о плотности твердого и жидкого германия до 1900 К [4, 5, 9, 39], сжимаемости [40], об энергии германия до 2000 К [6], вязкости [39, 41], коэффициентах самодиффузии [10, 33, 42], свойствах при высоких давлениях [11, 43, 44], ударном сжатии [12, 45] и др.
Рис. 5. ПКФ жидкого германия: 1 – дифракционные данные при 1253 К [1], 2 – модель с потенциалом ЕАМ. Невязка Rg = 0.020.
Значительное число работ посвящено моделированию жидкого германия с потенциалами типа SW [14, 15, 33, 46], Терсофа [18], а также с потенциалами сильной связи [44]. Ряд потенциалов приведен в репозитории NIST [22]. Потенциал МЕАМ также применили в [18, 24, 47], потенциал машинного обучения – в [48], метод ab initio в [44, 49]. Среднее КЧ жидкого германия (8.63) выше, чем у кремния, так что вклад ковалентной связи у германия меньше. Следовательно, потенциал ЕАМ для жидкого германия может оказаться вполне подходящим.
Парный вклад в потенциал ЕАМ для германия был определен нами алгоритмом Шоммерса, аналогично случаю кремния. ПКФ жидкого Ge при 1253 К рассчитали по структурному фактору [1] с включением метода наименьших квадратов [50]. Модель жидкости размером 2048 атомов была построена за 154 итерации с итоговой невязкой Rg = 0.020. Как видно из рис. 5, реальная и модельная ПКФ хорошо согласуются. Парный вклад в потенциал ЕАМ жидкого Ge показан на рис. 6 и приведен в табл. 5. Он также похож на жесткосферный со ступенькой вниз. Параметры потенциала погружения приведены в табл. 6.
Рис. 6. Парный вклад в потенциал ЕАМ германия, 1253 К.
Таблица 5. Парный вклад в потенциал жидкого германия. Алгоритм Шоммерса при 1253 К. Парная корреляционная функция [1] с поправкой [50], невязка Rg = 0.020
r, Å | ϕ(r), эВ | r, Å | ϕ(r), эВ | r, Å | ϕ(r), эВ | r, Å | ϕ(r), эВ |
1.5 | 122.300000 | 3.35 | –0.046844 | 5.2 | 0.002804 | 7.05 | 0.000684 |
1.55 | 101.303000 | 3.4 | –0.047240 | 5.25 | 0.001824 | 7.1 | 0.000647 |
1.6 | 82.303100 | 3.45 | –0.046399 | 5.3 | 0.000588 | 7.15 | 0.000611 |
1.65 | 65.292800 | 3.5 | –0.044678 | 5.35 | –0.001230 | 7.2 | 0.000504 |
1.7 | 50.263000 | 3.55 | –0.043242 | 5.4 | –0.003176 | 7.25 | 0.000431 |
1.75 | 37.220100 | 3.6 | –0.042177 | 5.45 | –0.005154 | 7.3 | 0.000489 |
1.8 | 26.207900 | 3.65 | –0.042683 | 5.5 | –0.006862 | 7.35 | 0.000379 |
1.85 | 17.311800 | 3.7 | –0.044350 | 5.55 | –0.008321 | 7.4 | 0.000296 |
1.9 | 10.611500 | 3.75 | –0.046536 | 5.6 | –0.009686 | 7.45 | 0.000160 |
1.95 | 6.066300 | 3.8 | –0.048615 | 5.65 | –0.010724 | 7.5 | 0.000095 |
2.00 | 3.391330 | 3.85 | –0.050613 | 5.7 | –0.011583 | 7.55 | –0.000021 |
2.05 | 2.031660 | 3.9 | –0.051642 | 5.75 | –0.012469 | 7.6 | 0.000084 |
2.1 | 1.338040 | 3.95 | –0.051700 | 5.8 | –0.013597 | 7.65 | 0.000363 |
2.15 | 0.847698 | 4.00 | –0.050500 | 5.85 | –0.014618 | 7.7 | 0.000651 |
2.2 | 0.407810 | 4.05 | –0.048100 | 5.9 | –0.015344 | 7.75 | 0.000941 |
2.25 | 0.242786 | 4.1 | –0.044900 | 5.95 | –0.015799 | 7.8 | 0.001186 |
2.3 | 0.163421 | 4.15 | –0.041400 | 6 | –0.015927 | 7.85 | 0.001263 |
2.35 | 0.098270 | 4.2 | –0.037700 | 6.05 | –0.015415 | 7.9 | 0.001323 |
2.4 | 0.045892 | 4.25 | –0.034200 | 6.1 | –0.014349 | 7.95 | 0.001494 |
2.45 | 0.002577 | 4.3 | –0.030800 | 6.15 | –0.013169 | 8 | 0.001467 |
2.5 | –0.033705 | 4.35 | –0.027300 | 6.2 | –0.011718 | 8.05 | 0.001491 |
2.55 | –0.063315 | 4.4 | –0.023400 | 6.25 | –0.009987 | 8.1 | 0.001490 |
2.6 | –0.086566 | 4.45 | –0.018900 | 6.3 | –0.008416 | 8.15 | 0.001489 |
2.65 | –0.103683 | 4.5 | –0.014400 | 6.35 | –0.007194 | 8.2 | 0.001543 |
2.7 | –0.114948 | 4.55 | –0.010100 | 6.4 | –0.006185 | 8.25 | 0.001735 |
2.75 | –0.120340 | 4.6 | –0.006120 | 6.45 | –0.005239 | 8.3 | 0.001765 |
2.8 | –0.120399 | 4.65 | –0.002540 | 6.5 | –0.004468 | 8.35 | 0.001859 |
2.85 | –0.115695 | 4.7 | –0.000011 | 6.55 | –0.003829 | 8.4 | 0.001948 |
2.9 | –0.107163 | 4.75 | 0.001530 | 6.6 | –0.003287 | 8.45 | 0.001727 |
2.95 | –0.095625 | 4.8 | 0.002330 | 6.65 | –0.002381 | 8.5 | 0.001419 |
3 | –0.082393 | 4.85 | 0.002780 | 6.7 | –0.001598 | 8.55 | 0.001148 |
3.05 | –0.069123 | 4.9 | 0.002900 | 6.75 | –0.000845 | 8.6 | 0.000767 |
3.1 | –0.057489 | 4.95 | 0.002840 | 6.8 | –0.000322 | 8.65 | 0.000519 |
3.15 | –0.049449 | 5 | 0.002930 | 6.85 | 0.000049 | 8.7 | 0.000388 |
3.2 | –0.045410 | 5.05 | 0.003130 | 6.9 | 0.000385 | 8.75 | 0.000199 |
3.25 | –0.044939 | 5.1 | 0.003420 | 6.95 | 0.000652 | 8.8 | 0.000043 |
3.3 | –0.045923 | 5.15 | 0.003320 | 7 | 0.000703 | 8.85 | 0.000000 |
Таблица 6. Жидкий германий. Параметры потенциала погружения F(r)
p1 | 0.57450 | n | 2048 |
p2 | 0.80 | dr, Å | 0.05 |
ρ1 | 0.875 | r0, Å | 9.86 |
ρ2 | 1.125 | rс, Å | 8.96 |
a1 | –2.7455, эВ | c1, эВ | 1.5559 |
При моделировании жидкого германия можно также использовать парный вклад в потенциал ЕАМ, заданный не в виде таблицы данных, а в виде кусочно-непрерывного полинома (7). В случае германия выбрали n = 7, k = 4. Коэффициенты в формуле (7) для германия приведены в табл. 7.
Таблица 7. Коэффициенты разложения парного вклада в потенциал ЕАМ германия
aim | Номер интервала i / Границы интервала ri – ri+1, Å | |||
1 / 2.10–2.80 | 2 / 2.80–4.20 | 3 / 4.20–6.00 | 4 / 6.00–8.95 | |
ai0 | –0.12914022803307D+00 | 0.31741657294333D-02 | –0.34244309645146D-02 | 0.00000000000000D+00 |
ai1 | 0.36918279528618D+00 | 0.12768036127090D+00 | 0.16339456196874D-02 | 0.00000000000000D+00 |
ai2 | –0.78049397398888D+01 | –0.12277981456606D+00 | 0.80042188493656D-01 | –0.29761810826526D-01 |
ai3 | –0.10718496961367D+03 | –0.65316546093368D+00 | 0.55860454460684D+00 | –0.12030635388921D+00 |
ai4 | –0.56677608973579D+03 | –0.15237788383222D+01 | 0.11235902756831D+01 | –0.15270893689047D+00 |
ai5 | –0.14389761345321D+04 | –0.20467986563964D+01 | 0.90001319360143D+00 | –0.86186740805348D-01 |
ai6 | –0.17372676812327D+04 | –0.12959413204484D+01 | 0.31329384464215D+00 | –0.22597590985478D-01 |
ai7 | –0.80443673950751D+03 | –0.29391029780645D+00 | 0.39797389267772D-01 | –0.22478041469157D-02 |
Этот потенциал работает в интервале 2.20–8.95 Å. К нему надо гладко добавить восходящую ветвь при r < 2.20 Å.
С помощью потенциала ЕАМ была построена серия моделей жидкого германия размером 2048 атомов в основном кубе с периодическими граничными условиями, при температурах 1200–1800 К, в режимах NVT и NpT. Парная корреляционная функция g(r) при 1253 К показана на рис. 5 в отличном согласии с дифракционными данными [1].
На рис. 7 показано распределение КЧ атомов жидкого германия при 1253 К, рассчитанное с радиусом сферы ближайших соседей 3.20 Å. Оно имеет вид, характерный для простых металлов, и проходит через максимум при Z = 6. Среднее значение Z = 6.17 ± 1.35. При выборе радиуса сферы соседей 3.60 Å получается Z = 8.63 ± 1.41. Как и в случае кремния, преимущественно тетраэдрическая упаковка не наблюдается. Распределение КЧ в германии очень похоже на случай кремния.
Рис. 7. Координационные числа модели германия при 1253 К. Радиус сферы ближайших соседей 3.6 Å.
Аналогично кремнию выглядит и распределение азимутальных углов (рис. 8), но со сдвигом максимума от величины 60° у кремния к значению 56°. При меньшем радиусе сферы 2.8 Å первый максимум находится при угле 65°. В методе ab initio [51] распределение углов аналогично, но при радиусе 2.8 Å первый максимум при 80° значительно ниже второго (при 100°). До тетраэдрического угла 109° все еще довольно далеко.
Рис. 8. Азимутальные углы в модели Ge при 1253 К: 1 – радиус сферы ближайших соседей 3.6 Å, 2 – радиус сферы 2.8 Å.
Поведение структуры жидких Si и Ge по сравнению с картиной при потенциалах SW и Ter указывает на ослабление ковалентной связи при смещении элементов сверху вниз по Периодической системе – от кремния к германию.
Результаты расчетов свойств жидкого германия в режиме NpT при температурах до 2000 К приведены в табл. 8. Литературные данные по плотности демонстрируют значительный разброс. Наилучшее согласие получено с данными [5, 52]. Видно также хорошее согласие с опытом по форме ПКФ при 1253 К (небольшая невязка Rg = 0.020), а также с ab initio ПКФ при 2000 К [51]. Расхождения с опытом по энергии моделей составляют всего 3 кДж/моль при 2000 К (без учета энергии электронов). Теплоемкость жидкого германия близка к 23 Дж/(моль К), то есть к классической величине 3R. Реальная теплоемкость жидкого германия немного выше – 27.6 Дж(моль К).
Таблица 8. Свойства моделей Ge, построенных методом МД при p ∼ 0.001 ГПа
T, K | d, г/см3 | <ρ>b | Rg | –E, кДж/моль | KT, ГПа | D×105, см2/с | ||||||
МД | Эксп [52] | Эксп [5] | Эксп [39] | –EMD | –Eexp [6] | МД | Эксп [31] | MD | Эксп [10] | |||
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 |
0 | – | – | – | – | – | – | – | 369.04 | – | – | – | – |
298a | – | – | 5.372 | – | – | – | 236.01 | 364.40 | – | – | – | – |
673 | – | – | 5.333 | – | – | – | – | – | – | – | – | – |
1153 | – | – | 5.284 | – | – | – | – | – | – | – | – | – |
1200 | 5.542 | – | – | 5.58 | – | – | 277.44 | 341.05 | – | – | 10.1 | 13.0 |
1210.4sol | – | – | – | – | – | – | 340.76 | – | – | – | – | |
1210.4liq | 5.526 | – | – | – | 1.00741 | – | 303.67 | 303.73 | – | – | – | – |
1253 | 5.510 | 5.47 | 5.49 | 5.556 | 1.00467 | 0.019 | 302.75 | 302.72 | 30.8 | 31.5 | 11.9 | 13.9 |
1300 | 5.485 | 5.45 | 5.46 | 5.535 | 1.00070 | 0.021 | 301.54 | 301.26 | – | 31.3 | 12.5 | 14.6 |
1400 | 5.436 | 5.41 | 5.43 | 5.487 | 0.99540 | – | 299.08 | 298.50 | – | 30.4 | 14.0 | 16.2 |
1500 | 5.387 | 5.35 | – | 5.445 | 0.98453 | – | 296.82 | 295.74 | 28.2 | 29.7 | 15.8 | 17.6 |
1600 | 5.342 | 5.30 | – | 5.400 | 0.97576 | – | 294.46 | 292.98 | – | – | 17.3 | 19.0 |
1700 | 5.295 | 5.26 | – | 5.352 | 0.96786 | – | 292.17 | 290.02 | – | – | 18.6 | 22 [51] |
1800 | 5.236 | 5.22 | – | 5.32 | 0.95652 | – | 289.72 | 287.49 | – | – | 19.5 | – |
1900 | 5.178 | 5.16 | – | – | 0.94962 | – | 287.37 | 284.70 | – | – | – | – |
2000 | 5.135 | – | – | – | 0.94245 | – | 285.18 | 281.94 | 18.7 | – | – | 23 [51] |
a Кристалл.
b Стандартные отклонения растут сверху вниз от 0.068 до 0.072.
Скорость звука в жидком германии измерена в [31] при температурах 1215–1443 К. При 1253 К адиабатическая сжимаемость равна βs = 2.52×10–11 м2/Н, откуда адиабатический модуль сжатия Ks = 1/βs = 39.7 ГПа. По данным МД, отношение теплоемкостей Cp/Cv = 1.26. Соответственно, изотермический модуль равен КТ = KsСv/Cp = 31.5 ГПа. При значении c1 = 1.5559 величина модуля модели, найденная по зависимости давления от объема, равна КТ = 30.8 ГПа. При нагревании от 1690 до 2000 К модуль КТ модели уменьшается до 18.7 ГПа.
Близость потенциала жидкого германия к жесткосферному (см. рис. 6) позволяет проверить формулу (5) для коэффициента самодиффузии. Принимая диаметр жесткой сферы σ = 2.45 Å, находим для модели при 1253 К N/V = 0.045719 ат/Å3, n = 0.6723, D0 = (3/8)×2.45×10–8×(8.31×107×1253/π/72.59)1/2 = 1.9632×10–4 см2/с, и, наконец, DHS = 1.13×10–4 см2/с. Эта величина хорошо согласуется с цифрой 1.19×10–4, полученной методом МД. Расчеты методом ab initio дают при 1250 К близкое значение D = 0.95×10–4 см2/с [49]. Псевдопотенциальный расчет дает 1.27×10–4 [53]. Результаты прямых измерений самодиффузии в жидком германии близки к приведенным выше (см. табл. 8).
В итоге получаем, что свойства жидкого германия хорошо описываются потенциалом ЕАМ, близким к жесткосферному со ступенькой вниз.
ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ
Ранее было показано методом ab initio, что при высоких давлениях химическая связь между двумя атомами в жидком кремнии возникает спонтанно вследствие случайного сближения этих атомов на расстояние меньше 2.5 Å [21]. В настоящей работе показано, что поведение моделей жидких кремния и германия хорошо описывается потенциалами ЕАМ с практически полным исчезновением направленности связи. При значениях координационных чисел 6 и выше более выгодным оказывается локально изотропный ближний порядок. В итоге гибридизация электронных состояний типа sp3 в жидких кремнии и германии не реализуется. С концепцией изотропного взаимодействия согласуются многие свойства жидких кремния и германия. Отсюда следует, что потенциалы, генерирующие направленность связи в кристаллах (SW, Ter и т. п.), должны быть мало пригодны для моделирования жидких кремния и германия при плотностях, близких к обычным. Этим объясняется относительно невысокая точность потенциалов, предлагаемых для описания одновременно твердых и жидких фаз кремния и германия.
Применим это рассуждение, например, для углерода, расположенного над кремнием в Периодической системе. При переходе вверх по 4-й группе элементов возрастают потенциалы ионизации атомов (от 7.90 у германия до 8.15 эВ у кремния и 11.3 у углерода) и, соответственно, увеличиваются границы температур и давлений, при которых должны происходить с ростом плотности постепенные переходы от направленной связи в углероде (типа sp, sp2 или sp3) к изотропной структуре жидкого углерода. Такой переход наблюдается при моделировании углерода методом ab initio при 9000 К и плотностях свыше 5.8 г/см3 [54]. Оставаясь в рамках классической молекулярной динамики, следует ожидать плавного изменения межчастичного потенциала при сжатии от варианта направленной связи (типа SW) к изотропному потенциалу. В области перехода его можно реализовать, например, в виде суммы этих двух потенциалов с весами, зависящими от плотности жидкости.
Ориентировочно, при понижении плотности в 2–3 раза (из-за нагревания) упомянутый переход может происходить в кремнии и германии в обратном направлении, а именно от изотропной жидкости вблизи от точки плавления – к жидкости с направленной связью. Это предсказание легко проверить методом ab initio.
Об авторах
Д. К. Белащенко
Национальный исследовательский технологический университет МИСИС
Автор, ответственный за переписку.
Email: dkbel75@gmail.com
Россия, Москва
Список литературы
- Waseda Y. The Structure of Non-Crystalline Materials. Liquids and Amorphous Solids. N.Y.: McGraw-Hill, 1980.
- Funamori N., Tsuji K. // Phys. Rev. Lett. 2002. V. 88. P. 255508.
- Demchuk T., Bryk T., Seitsonen A.P. et al. // arXiv:2009.00834. https://doi.org/10.48550/arXiv.2009.00834
- Assael M.J., Armyra I.J., Brillo J. et al. // J. Phys. Chem. Ref. Data. 2012. V. 41. № 3. https://doi.org/10.1063/1.4729873
- Глазов В.М., Чижевская С.Н, Глаголева Н.Н. Жидкие полупроводники. М.: Наука, 1967.
- Гурвич Л.А., Вейц И.В., Медведев В.А. и др. Термодинамические свойства индивидуальных веществ. Т. 2. Кн.2. М.: Наука, 1979.
- Desai P.D. // J. Phys. Chem. Ref. Data. 1986. V. 15. № 3. P. 967.
- Текучев В.В. Акустические и физико-химические свойства электронных расплавов. Волгоград. 2016.
- Регель А.Р., Глазов В.М. Физические свойства электронных расплавов. М.: Наука, 1980.
- Weis H., Kargl F., Kolbe M. et al. // J. Phys.: Condens. Matter. 2019. V. 31. P. 455101.
- Luo Sheng-Nian, Ahrens T.J., Asimow P.D. // J. Geophys. Res. 2003. V. 108. № B9. P. 2421. https://doi.org/10.1029/2002JB002317.
- Oleynik I.I., Zybin S.V., Elert M.L., White C.T. // CP845 “Shock Compression in Condensed Matter”. Ed. M.D. Furnish et al. 2005. P. 413.
- Stillinger F.H., Weber T.A. // Phys. Rev. B. 1985. V. 31. P. 5262.
- Dziedzic J., Principi E., Rybicki J. // J. Non-Cryst. Solids. 2006. V. 352. P. 4232.
- Jadhav P.P., Dongale T.D., Vhatkar R.S. // AIP Conference Proceedings 2162, 020038 (2019). https://doi.org/10.1063/1.5130248
- Tersoff J. // Phys. Rev. B. 1988. V. 38. P. 9902; 1989. V. 39. P. 5566.
- Ishimaru M., Yoshida K., Motooka T. // Phys. Rev. B. 1996. V. 53. № 11. P. 7176.
- Cook S.J., Clancy P. // Phys. Rev. B. 1993. V. 47. P. 7686.
- Bazant M.Z., Kaxiras E. // Phys. Rev. Lett. 1996. V. 77. P. 4370.
- Luo J., Zhou Ch., Cheng Y., Liu L. // J. Crystal Growth. 2020. P. 1. https://doi.org/10.1016/j.jcrysgro.2020.125785
- Štich I., Car R., Parrinello M. // Phys. Rev. B. 1991. V. 44. P. 4262.
- NIST. IPS Interatomic Potentials Repository: www.ctcms.nist.gov/potentials/refs.html
- Baskes M.I. // Phys. Rev. B. 1992. V. 46. P. 2727. https://doi.org/10.1103/PhysRevB.46.2727
- Baskes M.I., Nelson J.S., Wright A.F. // Phys. Rev. B. Condens. Matter. 1989. V. 40. № 9. P. 6085. https://doi.org/10.1103/physrevb.40.6085
- Starikov S.V., Lopanitsyna N.Yu., Smirnova D.E., Makarov S.V. // Computational Materials Science. 2018. V. 142. P. 303.
- Starikov S., Gordeev I., Lysogorskiy Yu. et al. // Computational Materials Science. 2020. V. 184. P. 109891. https://doi.org/10.1016/j.commatsci.2020.109891
- Daw M.S., Baskes M.I. // Phys. Rev. B. 1984. V. 29. № 12. P. 6443.
- Schommers W. // Phys. Rev. A. 1983. V. 28. P. 3599.
- Belashchenko David K. Liquid Metals. From Atomistic Potentials to Properties, Shock Compression, Earth’s Core and Nanoclusters. NOVA Science Publushers. NY.
- Zhu Z.G., Liu C.S. // Phys. Rev. B. 2000. V. 61. № 14. P. 9322.
- Hayashi M., Yamada H., Nabeshima N., Nagata K. // Int. J. Thermophysics. 2007. V. 28. № 1. P. 83. https://doi.org/10.1007/s10765-007-0151-9
- Chelikowsky J.R., Troullier N., Binggeli N. // Phys. Rev. B. 1994. V. 49. P. 114.
- Yu W., Wang Z.Q., Stroud D. // Phys. Rev. B. 1996. V. 54. № 19. P. 13946.
- Белащенко Д.К // Журн. физ. химии. 2019. T. 93. № 6. C. 877.
- Speedy R.J. // Mol. Physics. 1987. V. 62. № 2. P. 509.
- Бeлaщeнкo Д.K. // Physics–uspekhi. 2013. V. 183. № 12. P. 1176.
- Alteholz Th., Hoyer W. // J. Non-Cryst. Solids. 1999. V. 250–252. P. 48.
- Petkov V., Takeda S., Waseda Y., Sugiyama K. // J. Non-Cryst. Solids. 1994. V. 168. P. 97.
- Sato Y., Nishizuka T., Tachikawa T. et al. // High Temperatures – High Pressures. 2000. V. 32. P. 253.
- Tsuchiya Y. // J. Phys. Soc. Japan. 1991. V. 60. № 1. P. 227.
- Masaki T., Itami T. // “Modeling and Precise Experiments of Diffusion Phenomena in Melts under Microgravity” Annual Reports. 2002, NASDA-TMR-030005E.
- Kato M., Minowa S. // Trans. Iron Steel Institute of Japan. 1969. V. 9. P. 39.
- Tsuji K., Mori T., Hattori T. et al. // 2000B0087-CD-np BL04B1.
- Kōga J., Okumura H., Nishio K. et al. // Phys. Rev. B. 2002. V. 66. P. 064211.
- Kishimura H., Matsumoto H., Thadhani N.N. // J. Physics: Conference Series. 2010. V. 215. Р. 012145. https://doi.org/10.1088/1742-6596/215/1/012145
- Ding K., Andersen H.C. // Phys. Rev. B. 1986. V. 34. № 10. P. 6987.
- Kim Eun Ha, Shin Young-Han, Lee Byeong-Joo // Computer Coupling of Phase Diagrams and Thermochemistry. 2008. V. 32. P. 34.
- Zuo Y., Chen C., Li X. et al. // J. Phys. Chem. A. 2019. V. 124. № 4. P. 731. https://doi.org/10.1021/acs.jpca.9b08723
- Kresse G., Hafner J. // Phys. Rev. B. 1994. V. 49. № 20. P. 14251.
- Белащенко Д.К. // Кристаллография. 1998. Т. 43. № 3. С. 400.
- Kulkarni R.V., Aulbur W.G., Stroud D. // Phys. Rev. B. 1997. V. 55. P. 6896.
- Lucas L.D., Urbain G. // C. r. Acad. Sci. 1962. V. 255. № 19. P. 2414.
- Munejiri S., Shimojo F., Hoshino K., Itami T. // NASDA, Tsukuba 305–8505, Japan
- Hoshino K. // J. Phys.: Condens. Matter. 2009. V. 21. No 47. P. 474212. https://doi.org/10.1088/0953-8984/21/47/474212
Дополнительные файлы










