Журнал физической химии, 2023, T. 97, № 1, стр. 128-138

ВВЕДЕНИЕ

Свойства жидкого олова обсуждались во многих работах [1–5] в приближении парного взаимодействия. Однако значительный прогресс был достигнут при применении модели погруженного атома (Embedded Atom Model – ЕАМ [6]), на основе которой были построены модели большого числа кристаллических и ряда жидких металлов. Для олова также были предложены потенциалы ЕАМ и рассчитаны основные свойства моделей жидкости [7–9]. ЕАМ недостаточно хороша для описания температурных зависимостей свойств [10]. Недавно была предложена модифицированная модель ЕАМ–MEAM, в которую включен учет возможной угловой зависимости эффективной электронной плотности, создаваемой атомом в окружающем пространстве [11]. MEAM была применена, в частности, для расчета свойств Li, Ga, Sn, Be, Ni и др. В исходном варианте потенциала учитывались только ближайшие соседи данного атома. Далее было предложено учитывать не только ближайших соседей, но и следующих за ними – вариант 2NN MEAM (Ni, Sn, сплавы Sn–Pb [12–15], Li [16, 17] и др.).

В случае жидкого металла угловая зависимость потенциала маловероятна из-за изотропности жидкости. По данным [17], существенная разница расчетных свойств жидкого лития в вариантах 2NN MEAM и ЕАМ обнаружена не была. Поэтому ниже разработан более точный потенциал ЕАМ и дано сравнение результатов ЕАМ с МЕАМ [13] для жидкого олова.

Потенциалы EAM и MEAM. В модели ЕАМ потенциал модели погруженного атома имеет вид [6]:

Здесь U – потенциальная энергия, Φ(ρ_i) – потенциал погружения i-го атома, зависящий от эффективной электронной плотности ρ в месте нахождения центра атома, а вторая сумма по парам атомов содержит обычный парный потенциал φ(r). Эффективная электронная плотность в точке нахождения атома создается окружающими атомами и определяется по формуле ρ_i = $\sum\nolimits_j^{} {{{\psi }}({{r}_{{ij}}})} $, где ψ(r_ij) – сферически симметричный вклад в электронную плотность от соседа номер j. В расчетах используются три подгоночные функции Φ(ρ), φ(r) и ψ(r). Модель МЕАМ отличается от ЕАМ тем, что функция ψ(r) зависит от направления в пространстве и обеспечивает направленность химической связи в соответствующих (“плохих”) металлах [12].

Качество описания свойств жидких металлов на основе потенциалов ЕАМ и МЕАМ обсуждалось для случаев олова [13] и лития [17]. В [13] проведено моделирование олова в варианте 2NN MEAM и дано сравнение с результатами более ранней работы [12] с потенциалом того же типа. Ниже этот анализ проведен для олова с учетом более широкого набора свойств.

Парный вклад в потенциал ЕАМ. В настоящей работе потенциал ЕАМ для олова был рассчитан методикой, описанной в [9, 18, 19]. Вначале парный вклад в потенциал был найден с помощью алгоритма Шоммерса [4] по гистограмме дифракционной парной корреляционной функции (ПКФ) олова при 523 К [20]. Была применена, где это возможно, методика подавления ложных осцилляций ПКФ на малых расстояниях [21].

При сравнении графиков одного и того же свойства (например, двух различных гистограмм ПКФ g₁(r) и g₂(r)), будем определять степень их согласия как стандартное отклонение (“невязку”) с помощью формулы:

Здесь n₁ и n₂ – номера точек гистограммы ПКФ, между которыми вычисляется невязка. При невязке R_g < 0.04 две ПКФ визуально неразличимы. Так, в случае олова при 523 К невязка R_g между дифракционной [20] и модельной ПКФ равна всего 0.019.

Далее полученные гистограммы парных вкладов в потенциал Шоммерса аппроксимировали кусочно-непрерывными полиномами с 6 участками по оси расстояний (точки деления r₁, r₂, … r₇) по формуле:

Для олова выбрали k = 6, L = 8 и r_m = 2.80 Å. В этом выражении $H({{r}_{i}},{{r}_{{i + 1}}})$ – функция Хевисайда, равная 1 при r_i ≤ r ≤ ${{r}_{{i + 1}}}$ и нулю в остальных случаях, i – это номер интервала на оси r (i = 1, 2, …, 6). Условие непрерывности в точках r = r_i было применено к самому потенциалу и к его производной. Величина r₁ – это минимальное межчастичное расстояние в модели или немного меньшее значение. Радиус обрыва потенциала равен r_c = r₇ = 7.75 Å и является границей 2-й координационной сферы. В [12, 13] радиусы обрыва были меньше: 5.5 Å в [12] и 4.8 Å в [13]. Коэффициенты a_in парного вклада для расстояний 2.65 < r < 7.75 Å приведены в табл. 1. Они немного отличаются от предложенных в [8, 9].

При высоких давлениях атомы сближаются. Поэтому потенциал следует продолжить на область расстояний 0 < r ≤ 2.80 Å, где алгоритм Шоммерса при обычном давлении не работает из-за отсутствия нужных пар атомов. Вид потенциала на малых расстояниях можно определить по данным ударного сжатия. В итоге значения парного вклада на интервале 1.00–2.80 Ǻ рассчитываются интерполяцией данных алгоритма Шоммерса (табл. 2), а при r ≥ 2.80 Ǻ используется разложение (3). График парного вклада в потенциал ЕАМ показан на рис. 1.

Потенциал погружения ЕАМ. Потенциал погружения был аппроксимирован формулами [8, 9]:

причем ρ₀ = 1, а при ρ = ρ_i непрерывны сама функция Φ(ρ) и ее первая производная. Функция Φ(ρ) и все коэффициенты a, b и c выражаются в эВ. Координаты точек деления оси абсцисс возрастают в последовательности ρ₅–ρ₄–ρ₃–ρ₂–ρ₁–ρ₀–ρ₆–ρ₇–ρ₈. В итоге потенциал ЕАМ определяется параметрами p₁, p₂, a₁, c₁–c₉, ρ₁–ρ₈, m, n, q. Подгонку проводили по зависимости плотности и энергии олова от температуры вдоль бинодали, а также по данным статического и ударного сжатия. Выражения при ρ < ρ₀ используются при моделировании состояний с нормальной и пониженной плотностью, а при ρ > ρ₆ – для сжатых состояний. Параметр p₂ в (4) является подгоночным. Параметр p₁ определялся таким образом, чтобы получить для модели жидкости в “стандартном” состоянии (вблизи от точки плавления) среднее значение 〈ρ〉 = ρ₀ = 1. В этом случае потенциал погружения не влияет на движение частиц, поскольку dΦ(ρ)/dρ при ρ ≈ 1 очень мало. Коэффициенты a₂–a₉, b₂–b₉ рассчитываются из условия непрерывности потенциала погружения и его производной в точках ρ_i. Значения коэффициентов потенциала погружения приведены в табл. 3.

Проведение расчетов. Модели олова имели размер 2048 или 2000 атомов в основном кубе. Отдельные расчеты проводили на моделях большего размера. Моделирование проводили алгоритмом Л. Верле с шагом Δt = 0.01t₀, где единица времени t₀ = 1.109 × 10^–13 с. Применяли ансамбли NVT при реальной плотности олова [23] и NpT при давлении, близком к нулю. При расчетах учитывали электронные вклады в энергию и давление в модели свободных электронов (МСЭ), принимая четыре электрона на атом [7–9]. Значения тепловой энергии электронов E_eT приведены в табл. 4. Электронный вклад в давление p_eT рассчитывается по формуле p_eTV = (2/3)E_eT. Часть результатов опубликована ранее в [8, 9].

Результаты расчетов для бинодали олова. В табл. 5 и на рис. 2 показаны невязки R_g между модельными и фактическими ПКФ [20, 24]. При Т < 700 К невязки невелики (R_g ≤ 0.03), и обе ПКФ практически совпадают. При более высоких температурах расхождения между ПКФ увеличиваются (до R_g = 0.058 при 1573 К). Причинами расхождений являются как неточности потенциала ЕАМ, так и погрешность дифракционных данных, которая увеличивается с температурой. По данным графиков, в [13] невязка при 573 К довольно велика (R_g ≈ 0.12), хотя и несколько убывает при нагреве.

Результаты расчетов термодинамических свойств олова вдоль бинодали до 1973 К приведены в табл. 5. К энергии моделей были добавлены электронные вклады E_eT из табл. 4. Расчетная плотность отклоняется от фактической не более, чем на 1.9% (4 и 2 колонки). Электронные добавки в давление при невысоких температурах невелики и не превышают 0.1 ГПа. Согласие по энергии с данными опыта (10 и 11 колонки) в целом неплохое, но наблюдается небольшое занижение расчетных данных, которое постепенно увеличивается при нагреве. Это расхождение может быть отчасти следствием неточности МСЭ.

Далее был рассчитан с учетом электронных вкладов модуль всестороннего сжатия K_T, а также теплоемкости C_p и ${{C}_{{v}}}$ и скорость звука u по уравнению u = [(K_T/d)(C_p/${{C}_{{v}}}$)]^1/2, где d – плотность (см. табл. 6). Значения u приведены в табл.6 и показаны на рис. 3. Они согласуются с экспериментальными данными [25–27], а расхождения не превышают 3.6%.

Расчетный коэффициент самодиффузии монотонно увеличивается с температурой и описывается выражением D, см²/с = 5.5236 × 10^–9T^1.3344. Обзоры данных по самодиффузии олова приведены в [13, 22]. Характерен сильный разброс литературных данных. Значения D в условиях микрогравитации [28, 29] и близкие к ним значения, полученные в земных условиях [30], лежат заметно выше по отношению к нашим измерениям, значения в модели MEAM [13] также выше, а первопринципные значения [22] лежат немного ниже. Завышение значений D характерно для экспериментов с не полностью устраненной конвекцией.

Сдвиговая вязкость была рассчитана по Грину–Кубо через автокорреляционную функцию тензора вязких напряжений. С учетом работ [31, 32] вязкость η можно рассчитывать по уравнению:

Здесь k_B – постоянная Больцмана, и берется сумма по αβ = xy, xz, yz. Далее, где δ_αβ – символ Кронекера и Здесь m_j – масса атома j, r_jα и r_jβ – координаты атома  j, ${{{v}}_{{j\alpha }}}$ и ${{{v}}_{{j\beta }}}$ – компоненты его скорости, f_jkβ – β –компонента силы, с которой атом j действует на атом k. Расчеты проводили прогонами длиной 5000–10 000 шагов по времени, величины P_αβ(t) вычисляли на каждом шаге по времени, интегралы (5) рассчитывали на интервалах длиной 2500 шагов по времени. Усреднение путем сдвига расчетного интервала вдоль таблицы данных одного МД-прогона не проводилось. Расчетные значения вязкости при температурах 523–1173 К приведены в табл. 7. При всех температурах наблюдается занижение расчетных данных по отношению к опыту [23] на 10–30%, которое увеличивается с понижением температуры. Значения вязкости по Грину–Кубо, полученные с потенциалом МЕАМ [13], наоборот, завышены по отношению к данным [23] на 8–20%.

Рассмотрим, как выполняется соотношение Стокса–Эйнштейна, связывающее коэффициент самодиффузии и вязкость: D = kT/4πηr_a . Здесь r_a – “радиус атома”. При Т = 523 К по данным табл. 7 получаем r_a = 2.10 Ǻ. При учете данных МЕАМ [13] из табл. 7 получается r_a = 1.81 Ǻ. Если принять за диаметр атома координату 1-го пика ПКФ (3.17 Ǻ), то обе приведенные оценки завышены на 15–30%.

Состояния олова при высоких давлениях. Для расчета свойств олова при высоких давлениях требуется определить коэффициенты ρ₆, ρ₇, ρ₈, a₇–a₉, b₇–b₉, c₇–c₉, m, n, q. Их можно найти двумя способами: по форме ударной адиабаты и по данным статического сжатия при 298 К. При этом параметры потенциала ЕАМ, найденные при расчетах бинодали, сохраняются в обоих случаях.

Для расчета свойств в условиях ударного сжатия были использованы экспериментальные данные [33–35]. Ударная адиабата олова показана на рис. 4. Принято обозначение Z = V₀/V, где стандартный объем олова V₀ = 16.262 см³/моль. При давлениях свыше 310 ГПа имеется мало данных.

При высоких давлениях расстояния между ближайшими атомами убывают, так что необходима корректировка парного вклада в потенциал ЕАМ на малых расстояниях. Для уточнения парного вклада оказались полезными данные по ударному сжатию. В работе [38] было исследовано ударное сжатие олова и рассчитаны дифракционные ПКФ жидкого олова в нескольких состояниях на ударной адиабате (при давлениях 0, 52, 79 и 84 ГПа). Кроме того, в [38] провели моделирование олова методом ab initio в состояниях на ударной адиабате при давлениях 50.5, 76.3 и 91.8 ГПа. В [38] приведены графики ПКФ жидкого олова в условиях ударного сжатия.

Две важных характеристики формы ПКФ – минимальное межчастичное расстояние и высота 1-го пика ПКФ – позволяют подобрать крутизну парного вклада в потенциал ЕАМ на малых расстояниях. Именно этот парный вклад в потенциал ЕАМ приведен на рис. 1 и в табл. 2. Далее можно по форме ударной адиабаты найти описанным в [39] методом коэффициенты потенциала погружения, ответственные за поведение металла при высоких давлениях. Они приведены в табл. 3. При расчетах используется уравнение ударной адиабаты [40]:

Здесь E – энергия металла в сжатом состоянии, E₀ – энергия в исходном состоянии, p и V – давление и объем в сжатом состоянии, а p₀ и V₀ – они же в исходном состоянии. Величина E₀ = = –295.0 кДж/моль [7, 9]. Расчетные значения давления и энергии на адиабате с учетом электронных вкладов (в МСЭ при 4 эл/атом) показаны на рис. 4 и в табл. 8. Видно очень хорошее согласие с опытом для давления (табл. 8, колонки 2 и 11) и энергии (колонки 8 и 9) жидкого олова на адиабате. Электронные вклады быстро возрастают с увеличением температуры: до 238.59 кДж/моль и 19.56 ГПа при Z = 2. Температура на адиабате всюду ниже, чем в предположении, что электроны не дают вклада в энергию [41].

Дополнительное свидетельство адекватности потенциала ЕАМ можно получить, сравнивая ПКФ наших моделей с найденными методом ab initio и с ПКФ реального олова на ударной адиабате при Z = 1.5302 и 1.580 [38]. Эти ПКФ хорошо согласуются между собой. Высоты первых пиков дифракционных ПКФ в [38] несколько занижены (на ~0.5) по отношению к расчетам ab initio.

Холодное давление олова. Можно сравнить график холодного давления (изотерма при 298 К) с потенциалом ЕАМ с данными статического сжатия ОЦК олова при давлениях до ~200 ГПа [43–45] (рис. 5). При давлениях до 100 ГПа здесь получено очень хорошее согласие с опытом, которое довольно редко бывает при обработке данных ударного сжатия металлов. Однако при сжатии олова в 2 раза расхождение МД-расчета с опытом по давлению достигает 30 ГПа.

Нанокластеры олова. Поверхностные свойства нанокластеров анализировали ранее методом МД в ряде работ ([46–48] и др.). В частности, проводилась проверка применимости к нанокластерам макроскопических уравнений термодинамики (уравнения Лапласа для давления, формулы Толмена для зависимости поверхностного натяжения от кривизны поверхности, уравнения Кельвина для давления пара). В [46] исследовали методом МД нанокластеры с потенциалом Леннард-Джонса, а в [47] – нанокластеры серебра, построенные с потенциалом ЕАМ [5], а также методом ab initio. В [48] был исследованы нанокластеры Ar, а также Ag, Fe и Zn с потенциалами ЕАМ.

Поверхностное натяжение на границе жидкость–газ рассчитывают обычно, используя двухфазные МД-модели с плоской границей раздела. Погрешность такого расчета составляет ~15–25% (Li [17], Sn [13]), причем в сторону как завышения, так и занижения. Значительно проще определять методом МД не поверхностное натяжение σ, а избыточную энергию поверхности h. В работе [48] была предложена схема расчета поверхностной энергии сферических нанокластеров. В качестве объектов были использованы икосаэдрические нанокластеры Маккея [49], содержащие от 13 до 5233 атомов. Было показано, что зависимость энергии нанокластера E от числа атомов в нем N может быть с высокой точностью аппроксимирована уравнениями:

Коэффициенты a и b могут быть найдены графически. Значения поверхностной энергии E_s = h = = bN^2/3 были получены ранее для Ar, Ag, Fe, Zn [48], In [50] и Tl [51].

В силу термодинамического соотношения

можно было бы ожидать, что σ < h (поскольку обычно dσ/dT < 0), причем разница между σ и h должна быть невелика. В [47] величины σ и h оценивали для моделей кристаллических кластеров Ag размером от 13 до 5233 атомов. При расчетах методом ab initio (SIESTA) отклонения σ от h были невелики (разница в несколько процентов) и знакопеременны, а при расчетах методом МД с потенциалом ЕАМ из [5] величина σ при всех размерах кластеров Ag от 249 до 5233 атомов была на несколько процентов больше, чем h [47].

Аналогично [48], можно построить серию икосаэдрических нанокластеров Маккея для олова. При расчетах с потенциалом ЕАМ принималось во внимание, сколько атомов нанокластера находится не в жидкой, а в газовой фазе и не имеет ближайших соседей. Эти атомы при анализе кластеров не учитывались. Величины, рассчитанные без учета атомов в газовой фазе, обозначены как N*, E*, $E_{s}^{*}$, $S_{s}^{*}$. Модели находились в центральной части основного куба с длиной ребра 90 Å и не взаимодействовали со своими образами в соседних кубах. Периодически проводили остановку вращения кластера. Моделирование проводили прогонами по 10 000 шагов по времени. Шаг по времени равнялся Δt = 0.01t₀, где единица времени t₀ = 1.109 × 10^–13 с. Значения энергии этих нанокластеров были рассчитаны методом МД при 500 К и приведены в табл. 9 и на рис. 6. В последней строке приведены данные для модели сплошной жидкой фазы с периодическими граничными условиями (ПГУ).

Зависимость (7) для нанокластеров олова хорошо выполняется при N > 55, причем a = = ‒2.9146 эВ/атом, а b = 1.3044 эВ/атом^2/3 (см. рис. 6). Величина a почти совпадает с удельной энергией жидкого олова при 500 К (‒2.9212 эВ/атом). Коэффициент b отвечает за поверхностную энергию кластеров [48, 50]. В наших обозначениях поверхностная энергия кластера равна $E_{s}^{*}$ = b(N*)^2/3. Эти величины для кластеров с размерами 13 ≤ N ≤ 5083 приведены в табл. 9 (6 колонка).

Площадь поверхности нанокластера определяли путем разложения кластера на симплексы Делоне [48, 50]. Те грани кластера, которые с внешней стороны грани не имеют соседнего симплекса, являются поверхностными. Сумму площадей этих граней S_s можно отождествить с площадью поверхности кластера. Однако поверхность нанокластера напоминает при этом черепичную крышу. Поэтому была введена поправка на сглаживание поверхности, которая учитывает наклон грани по отношению к линии, соединяющей грань с центром кластера. Для учета этой негладкости каждая поверхностная грань симплексов нормализовалась, т.е. поворачивалась перпендикулярно к вектору, соединяющему эту грань с центром масс кластера. Новая сумма нормализованных площадей граней симплекса (повернутых) $S_{s}^{*}$ оказывается ниже исходной в среднем на 3–20% в зависимости от числа атомов в кластере и его рыхлости. Удельная поверхностная энергия увеличивается при этом на соответствующую долю.

В табл. 9 приведены значения нормализованной поверхности кластеров $S_{s}^{*}$ и удельной поверхностной энергии кластеров h = $E_{s}^{*}$/$S_{s}^{*}$. На рис. 7 показана зависимость удельной поверхностной энергии h = $E_{s}^{*}$/$S_{s}^{*}$ от размеров нанокластера. Здесь также величина h линейно зависит от N^–1/3 при N ≥ 55. В случае олова при 500 К предельное значение h = $E_{s}^{*}$/$S_{s}^{*}$ при N → ∞ (то есть в макроскопическом пределе) равно 0.416 Дж/м².

Вернемся к соотношению (8). В случае реального олова при 500 К σ = 0.552 Дж/м² и (dσ/dT)_p = = –0.13 мДж/(м² K) [52], так что при 500 К получаем h = 0.617 Дж/м². Эта величина выше полученного МД-расчетом значения 0.416 Дж/м². В отличие от условия σ < h, здесь получается обратное неравенство σ > h.

В работе [13] авторы измеряли методом МД с потенциалом MEAM поверхностное натяжение олова на плоской границе “модель олова–вакуум” и получили при 500 К значение σ ~ 0.645 Дж/м², которое больше фактического (0.552) и больше, чем h.

Такой же анализ был проведен на моделях лития при 500 К с использованием потенциала ЕАМ из работы [53]. На серии нанокластеров Маккея получено нормализованное значение h = 0.620 Дж/м² (при использовании не нормализованной поверхности кластеров получается 0.530 Дж/м²). С потенциалом [53, 54] в работе [17] измерено поверхностное натяжение лития на плоской границе жидкость–пар двухфазным методом и получено σ ≈ 0.480 Дж/м². Здесь неравенство σ < h выполняется. Расчет по опытным данным [52] при 500 К дает σ = 0.390 Дж/м², и при (dσ/dT)_p = –0.18 мДж/(м² К) получается h = = 0.480 Дж/м². Таким образом, расчеты на плоской границе раздела и на сферических нанокластерах согласуются между собой, но расходятся с опытом [52].

В случае индия реальное значение поверхностного натяжения при 433 К известно с разбросом от 0.556 до 0.595 Дж/м² (в среднем 0.571), а производная (dσ/dT)_p = –0.11 мДж/(м² К) [52]. Для величины h при 433 К по уравнению (8) получается h = 0.619 Дж/м². В то же время рассчитанная на сферических нанокластерах удельная поверхностная энергия индия E_s/S_s в пределе N → ∞ равна 0.474 Дж/м² [50] и ниже реальной величины h на ~25%.

В случае жидкого таллия при 588 К поверхностное натяжение равно 0.450 Дж/м² и (dσ/dT)_p = = –0.119 мДж/(м² К) [52]. Отсюда по формуле (8) находим h = 0.520 Дж/м². В работе [51] на серии кластеров Маккея описанным выше методом получено значение h = 0.401 Дж/м², то есть на 23% ниже. Аналогичное занижение поверхностной энергии по отношению к величине σ было получено для кластеров Ag в [47, 48].

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Сравнение результатов ЕАМ и МЕАМ показывает, что переход от варианта ЕАМ к МЕАМ не приводит к улучшению согласия с экспериментом для жидкого олова. Направленность химической связи в жидком олове практически отсутствует, и нет необходимости прибегать к использованию МЕАМ.

Некоторые результаты МД-расчетов поверхностной энергии и поверхностного натяжения кластеров в пределе N → ∞ приведены в табл. 10. Эти расчеты приводят к значительным расхождениям не только с опытом, но и с уравнением (8). Вряд ли эти расхождения можно приписать только недостаткам потенциалов ЕАМ/МЕАМ. Видимо, сама схема ЕАМ не подходит для расчетов поверхностных свойств металлов. При последовательном рассмотрении поверхностных свойств методами теории металлов необходимо учитывать “двухкомпонентность” металла (ионы + + электроны) [55], которая приводит к “выпячиванию” в вакуум электронного заряда (spill out) над кристаллической решеткой металла и к образованию двойного электрического слоя на поверхности. Эта двухкомпонентность не сочетается с идеологией ЕАМ. Поэтому нет оснований рассчитывать на предсказательную силу ЕАМ/МЕАМ в отношении поверхностных свойств. Например, довольно прямолинейный учет двойного слоя в случае серебра сильно меняет значения σ и h [47].

Пока неясны перспективы проверки соотношения (8) для нанокластеров. Теория функционала плотности (DFT) позволяет рассчитывать удельную поверхностную энергию с ошибкой не менее 10–20% [55]. Такова же и ошибка расчета поверхностной энергии h и поверхностного натяжения σ методом молекулярной динамики или прямого эксперимента. Однако если МД-расчет поверхностной энергии имеет ясный физический смысл, то МД-расчет поверхностного натяжения на двухфазных моделях металла через компоненты тензора напряжений приводит к результатам, такого смысла не имеющим.