Журнал физической химии, 2022, T. 96, № 7, стр. 952-958
О применении теории жидкостей к системам с потенциалом модели погруженного атома
a Национальный исследовательский технологический университет “Московский институт стали и сплавов”
Москва, Россия
* E-mail: dkbel75@gmail.com
Поступила в редакцию 22.12.2021
После доработки 22.12.2021
Принята к публикации 10.01.2022
- EDN: YDQNIK
- DOI: 10.31857/S004445372207007X
Аннотация
Рассмотрен вопрос о выполнимости следствий из уравнения Боголюбова–Борна–Грина (ББГ) для жидкостей и аморфных тел с гибридными потенциалами определенного вида (“теоремы B”). На примере жидкого таллия показано, что в случае жидкости с многочастичным потенциалом модели погруженного атома (ЕАМ) теорема B выполняется при не слишком больших значениях параметра гибридизации, как и должно быть при парных потенциалах. Причина выполнимости заключается в том, что у жидкостей с плотной структурой потенциал ЕАМ при данной температуре эквивалентен некоторому эффективному парному потенциалу, подотчетному уравнению ББГ.
Как известно, классическая теория жидкостей [1–6] описывает системы с самыми разными потенциалами межчастичного взаимодействия. Однако большинство результатов теории относится к системам с парным взаимодействием, когда потенциал имеет вид φ(r), где r – межчастичное расстояние. Часто применяемый в настоящее время потенциал модели погруженного атома (Embedded atom model – EAM) не является парным, так как зависит от координат многих окружающих атомов. В этой модели потенциальная энергия системы описывается уравнением [7]:
Здесь первая сумма содержит Φ(ρi) – потенциал погружения i-го атома, зависящий от эффективной электронной плотности ρ в месте нахождения центра атома, а вторая сумма по парам атомов содержит обычный парный потенциал. “Эффективная электронная плотность” ρi в точке нахождения атома создается окружающими атомами и определяется по формуле: где ψ(rij) – вклад в эффективную электронную плотность атома i от соседа номер j. В целом ЕАМ – всего лишь вычислительная схема, в которой используются функции, не обязательно точно соответствующие своему названию.В расчетах используются три подгоночные функции Φ(ρ), φ(r) и ψ(r), так что возможности согласования расчетных свойств с экспериментальными значениями очень широки. В случае кристаллов, где набор ближайших межатомных расстояний в состоянии равновесия невелик, удается правильно подогнать к опытным данным плотность, энергию, упругие постоянные, энергию образования вакансии, поверхностные свойства и т.д., а также относительную устойчивость различных кристаллографических модификаций данного металла. Схема ЕАМ позволяет обобщить ее на случай двойных металлических систем. Потенциалы ЕАМ считаются “трансферабельными”, т.е. для пар 11, 22 и т.д. их можно переносить (по определенным правилам) на кристаллический раствор или соединение. Требуется лишь определить дополнительно парные потенциалы для перекрестных пар (12, 13, 23…).
Если потенциал погружения не является просто числом, то энергия U не равна просто сумме по парам атомов, и взаимодействие не будет парным. Поэтому ряд уравнений теории жидкостей, полученных в приближении парного взаимодействия, может оказаться некорректным в случае потенциалов ЕАМ. В частности, сюда относится уравнение, связывающее основную структурную функцию жидкости – парную корреляционную функцию (ПКФ) g(r) – с парным межчастичным потенциалом φ(r). Это уравнение Боголюбова–Борна–Грина (ББГ) имеет вид:
(3)
$\begin{gathered} kT{{\nabla }_{1}}\ln {{R}_{2}}({{{\mathbf{r}}}_{1}},{{{\mathbf{r}}}_{2}}) = --{{\nabla }_{1}}u(|{\kern 1pt} {{{\mathbf{r}}}_{1}}--{{{\mathbf{r}}}_{2}}{\kern 1pt} |)-- \hfill \\ - \;\int {{{\nabla }_{1}}u(|{\kern 1pt} {{{\mathbf{r}}}_{1}}--{{{\mathbf{r}}}_{3}}{\kern 1pt} |)} \frac{{R{\text{3(}}{\mathbf{r}}{\text{1,}}\,{\mathbf{r}}{\text{2,}}\,{\mathbf{r}}{\text{3)}}}}{{R{\text{2(}}{\mathbf{r}}{\text{1,}}\,{\mathbf{r}}{\text{2)}}}}d{\mathbf{r}}{\text{3}}{\text{.}} \hfill \\ \end{gathered} $(4)
${{R}_{3}}({{{\mathbf{r}}}_{1}},{{{\mathbf{r}}}_{2}},{{{\mathbf{r}}}_{3}})\sim {{R}_{2}}({{{\mathbf{r}}}_{1}},{{{\mathbf{r}}}_{2}}){{R}_{2}}({{{\mathbf{r}}}_{1}},{{{\mathbf{r}}}_{3}}){{R}_{2}}({{{\mathbf{r}}}_{2}},{{{\mathbf{r}}}_{3}}).$Как показано в работе [8], уравнение ББГ остается справедливым и при низких температурах вплоть до абсолютного нуля. В этом случае оно имеет смысл уравнения механического равновесия системы под действием внутренних сил (см. также [9]) и принимает вид:
(5)
$\begin{gathered} 0 = --{{\nabla }_{1}}u(|{\kern 1pt} {{{\mathbf{r}}}_{1}}--{{{\mathbf{r}}}_{2}}{\kern 1pt} |)-- \\ - \;\int {{{\nabla }_{1}}u(|{\kern 1pt} {{{\mathbf{r}}}_{1}}--{{{\mathbf{r}}}_{3}}{\kern 1pt} |)} \frac{{R{\text{3(}}{\mathbf{r}}{\text{1}},{\mathbf{r}}{\text{2}},{\mathbf{r}}{\text{3)}}}}{{R{\text{2(}}{\mathbf{r}}{\text{1}},{\mathbf{r}}{\text{2)}}}}d{\mathbf{r}}{\text{3}}{\text{.}} \\ \end{gathered} $В случае реализации парного потенциала справедливо следующее. Предположим, что теорема А выполняется и что существуют два парных межчастичных потенциала u1(r) и u2(r), первый из которых генерирует данную структуру при некоторой температуре Т > 0, а второй генерирует при Т = 0 ту же структуру с такой же плотностью и с теми же функциями g(r) и R3(r1, r2, r3). Тогда первый удовлетворяет уравнению ББГ (3) при Т > 0, а второй – уравнению (5) с нулевой левой частью. Тогда, очевидно, суммарный потенциал u(r) = = u1(r) + λu2(r) будет также удовлетворять уравнению ББГ с той же ПКФ при любом коэффициенте λ. Для устойчивости системы требуется лишь, чтобы величина λ была положительной.
Таким образом, должно существовать семейство гибридных парных потенциалов
которое может генерировать при различных λ одну и ту же структуру некристаллической (жидкой или аморфной) системы [10, 11]. Назовем этот результат теоремой B.Следует учесть, что теорема А может не выполняться при значительных преобразованиях структуры вследствие изменений температуры или плотности. Кроме того, ПКФ определяются из дифракционных экспериментов с некоторой неточностью. Все же можно полагать, что для простых жидкостей практически полное совпадение ПКФ при одинаковой плотности числа частиц обеспечивает хорошее согласие трехчастичных корреляционных функций. Тогда указанная выше теорема B будет верна, по крайней мере, при не очень больших λ.
За степенью согласия между двумя ПКФ моделей H1 и H2 можно следить, вычисляя стандартное отклонение (“невязку”) между ними по формуле:
(7)
${{R}_{g}} = {{\left\{ {\frac{{\text{1}}}{{n{\text{2}} - n{\text{1}} + {\text{1}}}}\sum\limits_{n{\text{1}}}^{n{\text{2}}} {{{{[{{g}_{2}}({{r}_{j}})--{{g}_{1}}({{r}_{j}})]}}^{2}}} } \right\}}^{{1/2}}}{\kern 1pt} ,$Приведенные соображения поддаются экспериментальной проверке в компьютерном эксперименте. В [10] теорема B была проверена на примере системы с потенциалом
и на примере двухкомпонентной системы с парными межчастичными потенциалами Леннард-Джонса. Поскольку до 2016 г. не были предложены алгоритмы построения модели аморфной фазы по заданной ПКФ, то в [10] был использован следующий алгоритм. Вначале при Т = 0 построили модель аморфной фазы (модель H2) с потенциалом (8) методом непрерывной статической релаксации (см. ниже), а затем с помощью алгоритма Шоммерса [12] построили методом молекулярной динамики (МД) при 300 К модель H1 с такой же ПКФ, как у модели H2. В обозначениях теоремы B межчастичный потенциал модели H1 представляет собой u1 из (6), а потенциал модели H2 – соответственно u2. Далее была построена серия моделей (H3) при 300 К с потенциалами вида (6) и коэффициентами 1 < λ < 500. Согласно теореме B, ПКФ серии моделей H3 функция g(r) должна совпадать с ПКФ g1(r) моделей H1 и H2. За степенью согласия между двумя ПКФ моделей H1 и H3 следили, вычисляя стандартное отклонение (“невязку”) Rg между ними. Если величина Rg составляет не более нескольких сотых, то графики функций g(r) и g1(r) визуально практически неразличимы. В работе [10] было показано методом молекулярной динамики, что указанная выше теорема B о семействе гибридных потенциалов справедлива в пределах ошибки расчета. Например, при температуре 300 К и λ = 100 невязка Rg составляла всего 0.05–0.06. Это означает также (что очень важно), что теорема А о совпадении трехчастичных функций при одинаковых ПКФ справедлива даже вблизи абсолютного нуля.В настоящей работе проведена дополнительная проверка теоремы B методом МД на примерах жидкости с парным потенциалом, а также с потенциалом ЕАМ с учетом того, что в 2016 г. был предложен универсальный алгоритм построения модели аморфной фазы по известной ПКФ при Т = 0 [13]. В этом случае модель М1 строится при T > 0 с известным потенциалом φ1(r), а модель М2 строится при Т = 0 по ПКФ уже построенной модели М1. Следовало выяснить, выполняется ли теорема B в случае жидкостей с многочастичным взаимодействием – с потенциалом ЕАМ.
ПРОВЕРКА ТЕОРЕМ НА ПРИМЕРЕ ЖИДКОГО ТАЛЛИЯ
Исходная модель жидкого таллия при 1073 К. В работе [14] была построена серия моделей таллия в состояниях на бинодали при температурах 588–3000 К методом молекулярной динамики (МД) [15]. Потенциал ЕАМ этой серии (ниже обозначен как P0) был подобран по ПКФ жидкого таллия при 588 К и по зависимости плотности и энергии жидкого таллия от температуры. Для дальнейших расчетов выбрали модель таллия при 1073 К (модель М0). На рис. 1 показаны ПКФ реального таллия при 1073 К [16] в сравнении с ПКФ модели М0. Ее генерирует потенциал ЕАМ P0. Невязка между двумя ПКФ невелика (Rg = = 0.0516). При 588 К невязка еще меньше: Rg = = 0.0358.
Для проверки теоремы B в случае парных потенциалов требуется построить модель М1 с парным межчастичным взаимодействием при 1073 К и с ПКФ, совпадающей с ПКФ модели М0 жидкого таллия с потенциалом ЕАМ P0. Назовем моделями пара-таллия молекулярно-динамические модели жидкого или аморфного вещества, полученные из модели жидкого таллия вариацией межчастичного потенциала и температуры, но имеющие такую же, как у него, парную корреляционную функцию. Модель пара-таллия М1 c ПКФ реального жидкого таллия при 1073 К была построена при 1073 К алгоритмом Шоммерса. Этот алгоритм восстанавливает парный потенциал по заданной ПКФ. При этом удалось достигнуть очень низкой невязки Rg = 0.0083 между ПКФ моделей М0 и М1. При такой невязке график ПКФ модели М1 совпадает с графиком ПКФ модели М0 (кривая 2 на рис. 1). Полученный парный потенциал P1 (т.е. φ(r), не ЕАМ!!) и силовая функция F(r) = –dφ(r)/dr показаны на рис. 2.
Кроме модели М1 пара-таллия, для проверки теоремы B необходимо также построить модель аморфной фазы при Т = 0 с парным потенциалом и с ПКФ, совпадающей с ПКФ моделей М0 и М1. Модели аморфной фазы при Т = 0 можно строить методом непрерывной статической релаксации (НСР) [17]. Этот метод похож на метод МД, но смещение атомов на каждом шаге МД проводится на одну и ту же величину Δr в направлении равнодействующей силы. Предлагавшиеся ранее алгоритмы (например, [18, 19]) не позволяли получить хорошее усреднение функции ПКФ.
Если просто применить потенциал P1 и построить модель при Т = 0 с плотностью модели М1 10.6827 г/см3 методом НСР, то получается модель аморфной фазы (М2) с очень высоким первым пиком ПКФ (рис. 3). Характерная рябь на ПКФ обусловлена недостаточным усреднением в методе статической релаксации.
Искомую модель М3 при Т = 0 с ПКФ моделей М0–М1 можно построить с помощью универсального алгоритма сравнения координационных чисел [13]. Он заключается в том, что для исходной модели М1 рассчитываются средние числа соседей в координационных сферах с радиусами rm = = mΔr (m = 1, 2, 3, …, mc), где шаг Δr = 0.05 Å, 0 < < rm < rc, и радиус обрыва взаимодействия rc = = 7.65 Å. Для искомой модели М3 такой расчет проводится на каждой итерации алгоритма. Обозначим эти числа соседей через n1(m) и n2(m), а парную силу, действующую на атом i со стороны атома j, через F(rij). Числа n1(m) гистограммы координационных чисел от времени не зависят. На каждом шаге универсального алгоритма таблица значений парной силы F(rm) модели М3 корректируется по формуле:
(9)
$\Delta F({{r}_{{\text{m}}}}) = \alpha [{{n}_{2}}(m)--{{n}_{1}}(m)],\quad m = 1,2,3, \ldots ,{{m}_{{\text{c}}}}.$Коэффициент α подбирается так, чтобы добавки ΔF(rm) были невелики. В итоге парная сила F(rm) на заданном межчастичном расстоянии увеличивается, если в сфере радиуса rm число соседей модели М3 больше, чем в модели М1, и уменьшается в обратном случае. На каждой итерации алгоритма коррекция (9) выполняется на каждом расстоянии гистограммы rm.
Шаг смещений атомов указанного алгоритма равнялся 0.01 Å. Размер модели М3 составлял 2000 атомов в основном кубе. Плотность модели равнялась 10.6827 г/см3. За 162 итерации была построена методом НСР модель М3 аморфного пара-таллия при 0 К с невязкой между ее ПКФ и ПКФ жидкого таллия при 1073 К с той же плотностью (модели М0–М1) всего Rg = 0.0314 (рис. 4). Ввиду малости Rg графики ПКФ моделей М1 и М3 практически совпадают.
На рис. 5 показан парный потенциал φ(r) модели М3 пара-таллия при Т = 0, рассчитанный интегрированием силовой функции F(r) = –dφ(r)/dr. Назовем этот потенциал (парный) P3. Следует учитывать, что рассчитанные при Т = 0 силовая функция и потенциал определяются с точностью до умножения на любое положительное число [20].
ПРОВЕРКА СПРАВЕДЛИВОСТИ ТЕОРЕМЫ B В СЛУЧАЕ ПАРНЫХ ПОТЕНЦИАЛОВ
Итак, построены модель жидкого пара-таллия (М1) при 1073 К и плотности 10.6827 г/см3 с парным потенциалом P1 и модель аморфного пара-таллия (М3) при Т = 0 К с парным потенциалом P3. ПКФ этих моделей практически совпадают с ПКФ жидкого таллия при 1073 К. Для проверки теоремы B следует провести моделирование жидкостей с гибридными потенциалами вида:
где φ1 – парный потенциал модели жидкого пара-таллия М1, а φ3(r) – парный потенциал P3 аморфного пара-таллия (модель М3), восстановленный универсальным алгоритмом при 0 К и показанный на рис. 5. Расчеты проводили методом МД на моделях размером 2000 атомов в основном кубе при постоянной температуре 1073 К и плотности 10.6827 г/см3. Применяли алгоритм Л. Верле. Результаты приведены в табл. 1.Таблица 1.
λ | Rg | р, ГПа | –Е, кДж/моль | D × 105, см2/с |
---|---|---|---|---|
0 | 0.0000 | –0.0576 | 61.99 | 6.03 |
0.2 | 0.0054 | –0.0381 | 61.78 | 6.27 |
0.4 | 0.0051 | –0.0234 | 61.50 | 6.05 |
0.7 | 0.0054 | 0.0131 | 61.10 | 5.86 |
1.0 | 0.0054 | 0.0444 | 60.71 | 5.97 |
2.0 | 0.0074 | 0.1471 | 59.34 | 6.00 |
5.0 | 0.0164 | 0.4288 | 55.40 | 5.92 |
10.0 | 0.0286 | 0.9391 | 48.70 | 5.84 |
20.0 | 0.0546 | 1.9444 | 35.44 | 5.53 |
Из табл. 1 видно, что невязка между ПКФ модели М4 с гибридным потенциалом и ПКФ модели с λ = 0 (т.е. с ПКФ модели М1) остается очень низкой даже при λ = 10, хотя энергия при этом уже на 20% выше, чем при λ = 0. С ростом λ коэффициент самодиффузии очень слабо убывает. На рис. 6 показана ПКФ М4 модели с λ = 10.
Это означает, что теорема B выполняется вплоть до λ ~ 10. При более высоких значениях λ теоремы А и B уже не выполняются.
ПРОВЕРКА СПРАВЕДЛИВОСТИ ТЕОРЕМЫ B В СЛУЧАЕ ПОТЕНЦИАЛОВ ЕАМ
В этом случае гибридный потенциал может включать потенциал ЕАМ P0 модели жидкого таллия М0 при 1073 К и потенциал P3 модели М3. Построение моделей M5 с потенциалами
проводили методом МД при 1073 К и плотности 10.6827 г/см3. При этом в качестве φЕАМ выбрали потенциал ЕАМ модели таллия М0 из [14]. Потенциал φ3(r) выбрали так же, как в предыдущем случае (10). Модели содержали по 2000 атомов в основном кубе. При расчете невязок в качестве исходной выбрали ПКФ модели М0 жидкого таллия при 1073 К. Результаты приведены в табл. 2.Таблица 2.
λ | Rg | р, ГПа | –Е, кДж/моль | D × 105, см2/с |
---|---|---|---|---|
0 | 0.0047 | 0.0040 | 148.74 | 5.90 |
0.2 | 0.0059 | 0.0130 | 148.49 | 6.18 |
0.4 | 0.0069 | 0.0423 | 148.22 | 6.27 |
0.7 | 0.0050 | 0.0601 | 147.88 | 5.83 |
1.0 | 0.0060 | 0.0854 | 147.44 | 5.81 |
2.0 | 0.0076 | 0.1799 | 146.14 | 5.86 |
5.0 | 0.0137 | 0.4629 | 142.13 | 6.08 |
10.0 | 0.0256 | 0.9543 | 135.43 | 5.71 |
20.0 | 0.0495 | 1.9129 | 122.08 | 5.21 |
Из табл. 2 видно, что при увеличении λ вплоть до ~10 невязки Rg между ПКФ моделей М0 и М5 невелики, т.е. парная корреляционная функция таллия практически не изменяется (см. рис. 7), хотя давление и энергия заметно увеличиваются, а коэффициент самодиффузии убывает. Это означает, что выполняются теоремы А и B. При более высоких значениях λ теоремы А и B перестают выполняться.
Таким образом, при не слишком высоких значениях λ жидкости и аморфные системы с потенциалами ЕАМ ведут себя в отношении теоремы B аналогично системам с парными потенциалами.
ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ
Можно обсудить причины такого согласия между системами с парным и многочастичным характером межчастичного взаимодействия. Объединим уравнения (1) и (2) ЕАМ:
(12)
$U = \mathop \sum \limits_i \,\Phi \left[ {\mathop \sum \limits_{j \ne i} \,\psi ({{r}_{{ij}}})} \right] + \sum\limits_{i < j} {\varphi {\text{(}}rij{\text{)}}} .$Более заметные отклонения от теорем А и В следует ожидать в состояниях жидкости с более широким разбросом значений эффективной электронной плотности, т.е. при повышенных температурах и пониженных плотностях. Например, в случае модели таллия при 3000 К и близком к нулю давлении 〈ρ〉 = 0.5905 ± 0.1538 [14]. При таком широком разбросе величин ρi теорема В может перестать выполняться при меньших значениях λ.
Возвращаясь к универсальному алгоритму, следует иметь в виду, что задача построения модели аморфной фазы при Т = 0 с заданной ПКФ не всегда имеет решение. Например, в случае модели таллия при 3000 К не удается построить указанным алгоритмом аморфную модель при Т = 0 с ПКФ жидкого таллия при 3000 К. Плотность таллия при 3000 К и близком к нулю давлении равна 6.676 г/см3 и составляет всего 62.5% от плотности при 1073 К [14]. Поэтому, если бы такую модель удалось построить, то она была бы сильно растянута, имела бы отрицательное давление и оказалась бы неустойчивой относительно образования в ней внутренних полостей. Модель же с иной плотностью не отвечала бы условиям теорем А и В – идентичностью плотностей и структуры моделей с исходным и гибридным потенциалами.
Список литературы
Боголюбов Н.Н. Проблемы динамической теории в статистической физике. М.: Ростехиздат, 1946.
Физика простых жидкостей. Статистическая теория. Пер. с англ. Под ред. Г. Темперли, Дж. Роулинсона, Дж. Рашбрука. М.: Мир, 1971.
Крокстон К. Физика жидкого состояния. Статистическое введение. М.: Мир. 1978.
Hansen J.-P., McDonald I.R. Theory of Simple Liquids. Univ. Cambridge. Academic Press, 2006.
Куни Ф.М. Статистическая физика и термодинамика. М.: Наука, 1981.
Мартынов Г.А. Классическая статистическая механика. Теория жидкостей. Долгопрудный: Интеллект, 2011.
Daw M.S., Baskes M.I. // Phys. Rev. B. 1984. V. 29. № 12. P. 6443.
Белащенко Д.К., Менделев М.И. // Журн. физ. химии. 1995. Т. 69. № 3. С. 543.
Johnson M.D., Hutchinson P., March N.H. // Proc. Roy. Soc. (L). A. 1964. V. 282. P. 283.
Белащенко Д.К. // Журн. физ. химии. 2004. Т. 78. № 9. С. 1621.
Belashchenko D.K. Liquid Metals. From Atomistic Potentials to Properties, Shock Compression, Earth Core and Nanoclusters. Nova Science Publ., 2018.
Schommers W. // Phys. Rev. A. 1983. V. 28. P. 3599.
Белащенко Д.К. // Журн. физ. химии. 2016. Т. 90. № 4. С. 483.
Белащенко Д.К. // Там же. 2022. Т. 96. № 3. С. 390.
Norman G.E., Stegailov V.V. // Math. Models and Computer Simulations. 2013. V. 5. № 4. P. 305.
Waseda Y. The Structure of Non-Crystalline Materials. Liquids and Amorphous Solids. N.Y.: McGraw-Hill, 1980. 325 p.
Белащенко Д.К. // Физ. мет. и металловед. 1985. Т. 60. № 6. С. 1076.
Bennet Ch.H. // J. Appl. Phys. 1972. V. 43. C. 2727.
Boudreaux D.S., Gregor J.M. // J. Appl. Phys. 1977. V. 48. № 1. P. 152.
Белащенко Д.К. // Физ. мет. и металловед. 1987. Т. 63. № 4. С. 665.
Белащенко Д.К. // Журн. физ. химии. 2021. Т. 95. № 12. С. 1804.
Белащенко Д.К. // Там же. 2021. Т. 95. № 1. С. 80.
Дополнительные материалы отсутствуют.
Инструменты
Журнал физической химии