Журнал физической химии, 2019, T. 93, № 6, стр. 877-889

Сурьма принадлежит к небольшой группе аномальных металлов (полуметаллов), плавящихся с уменьшением объема. Данные о термохимических свойствах жидкой фазы (очень скудные, особенно в отношении теплоемкости и энтальпии) приведены в [1, 2], о плотности – в [3], о скорости звука и сжимаемости – в [4–6]. Коэффициенты самодиффузии жидкой сурьмы приведены в [7–10]. Структуру жидкой сурьмы исследовали неоднократно в [6, 11–19]. В дифракционных работах применяли как рентгеновское, так и нейтронное излучение. Характерный признак аномальной структуры сурьмы – наличие побочного максимума справа от первого пика парной корреляционной функции (ПКФ) (см. рис. 1). Таблицы ПКФ для двух температур приведены в [11, 12]. В [14] исследована структура жидкой сурьмы под давлением до 9.4 ГПа, при температурах на 50 К выше точки плавления.

Некоторые характеристики структуры жидкой сурьмы приведены в табл. 1. Виден значительный разброс координаты (r₁) первого пика ПКФ: от 2.83 до 3.28 Å. Эта координата является важной характеристикой структуры жидкости. У плотных некристаллических структур величина ρ₁ = = r₁(N/V)^1/3 = 1.08 ± 0.2 [20] (N – число частиц в объеме V). У сурьмы, с учетом данных [11, 12], при 923 К получается ρ₁ = 1.041; по данным [6], ρ₁ = = 0.967. У висмута, расположенного ниже в Периодической системе, при 553 К величина ρ₁ = 1.034 [11, 12] и 1.006 [21]. Таким образом, структура жидкой сурьмы (и в меньшей степени – висмута) рыхлая.

Работ по моделированию жидкой сурьмы сравнительно немного. В [24, 25] было проведено моделирование жидкой сурьмы методами Монте-Карло и молекулярной динамики (МД) с применением эффективных парных потенциалов. В [10, 23, 26] исследованы свойства жидкой сурьмы квантово-механическими методами. Результаты этих расчетов не совпадают. Так, в [23] высота первого пика ПКФ (g₁) сурьмы при 1072 К получена равной ∼2.3 и немного выше данных [11, 12] (2.17), а в [10] высота немного меньше 2.0, что согласуется с данными [6], но не с [23]. Кроме того, координаты r₁ в [10, 23] отличаются на 0.23 Å, а в [11, 16] – на 0.33 Å (табл. 1). Скорость звука в [10] превышает фактическую [5] в ∼1.24 раза, что означает завышение модуля всестороннего сжатия в ∼1.5 раза. Следовательно, силы межчастичного взаимодействия в [10] заметно более жесткие, чем в реальной сурьме.

Принято считать, что рыхлость структуры обусловлена участием направленной связи (как, например, в жидких Si и Ge), причем это участие ослабевает с ростом температуры. В рамках концепции парного взаимодействия это приводит к потенциалам, зависящим от температуры, что исключает возможность применения теории жидкостей и статистической механики вообще. Поэтому особый интерес представляет проверка возможности моделировать одноатомные жидкости с рыхлой структурой с помощью модели погруженного атома (Embedded Atom Model – EAM). Учет направленности связи в модифицированном варианте МЕАМ [27] в случае лития не привел к хорошим результатам [28].

При расчетах потенциала ЕАМ по методике [29] парный вклад в потенциал определяется методом МД с помощью алгоритма Шоммерса [30]. В этом алгоритме используется ПКФ жидкости, и правильный выбор целевой ПКФ очень важен. Вследствие различия значений r₁ эффективный диаметр частиц в [6] меньше, чем в [11] на 7%. Поскольку на начальном этапе наших расчетов не было оснований предпочесть какие-то из этих данных, то были использованы ПКФ работ [11] и [6] по отдельности. Для обоих вариантов можно получить межчастичные потенциалы, которые при 923–930 К приводят к хорошему согласию соответствующих расчетной и дифракционной ПКФ.

Потенциал ЕАМ. Мы применили потенциал ЕАМ в форме [31]:

где φ(r) – это обычный парный потенциал, а Φ(ρ) – потенциал погружения. Потенциал погружения выражается через эффективную электронную плотность ψ(r) и имеет вид [29, 32]:

Методика расчета потенциала ЕАМ для жидких металлов описана, например, в [20, 29, 32, 33]. Коэффициент a₁ определяется по энергии системы, c₁ – по модулю всестороннего сжатия, а остальные коэффициенты – по зависимости плотности от температуры вдоль бинодали и при высоких давлениях – по форме ударной адиабаты. Энергия рассчитывается по отношению к идеальному газу при абсолютном нуле, и для сурьмы в стандартном состоянии равна – 255.8 кДж/моль (так как для одноатомного газа сурьмы $\Delta H_{{298}}^{0}$ = 262 кДж/атом).

Расчет парного вклада в потенциал ЕАМ сурьмы. На первом этапе расчета применили алгоритм Шоммерса [30] с целевой ПКФ работы [6], полученной методом рассеяния нейтронов. При этом парный вклад в потенциал ЕАМ получается в виде гистограммы. На рис. 1 показаны дифракционная [6] и модельная ПКФ сурьмы при 923 К и реальной плотности 6.469 г/см³ [3] (давление модели близко к нулю: ∼10^{– 4} ГПа). Стандартное отклонение между этими кривыми (“невязка” R_g) составляет всего 0.035, и они визуально практически совпадают. Расчетная высота 1-го пика (∼2.00) хорошо согласуется с дифракционными данными. При 948 К согласие МД-расчетов с дифракционными данными также очень хорошее (рис. 1). Минимальное межчастичное расстояние при 923 К равно 2.50 Å [6]. Отметим, что поведение ПКФ на малых расстояниях определяется в дифракционных исследованиях довольно неточно в связи с особенностями фурье-преобразования структурного фактора. Полученную алгоритмом Шоммерса гистограмму парного вклада в потенциал аппроксимировали кусочно-непрерывным полиномом восьмой степени с семью участками по оси расстояний (точки деления r₁, r₂, …, r₈) по формуле:

Для сурьмы выбрали k = 7 и L = 8. В этом выражении H(r_i, r_i + 1) – функция Хевисайда, равная единице при r_i ≤ r ≤ r_{i + 1} и нулю в остальных случаях, i – это номер интервала на оси r (i = 1, 2, …, 7). Потенциал и его производная непрерывны в точках r = r_i. Параметры парного потенциала для расстояний 2.45 Å < r < 8.50 Å приведены в табл. 2. Радиус обрыва потенциала равен 8.50 Å. Этот потенциал был продолжен на область расстояний 0 < r ≤ r₁ = 2.45 Å (где алгоритм Шоммерса не работает) по формуле:

Здесь r выражено в Å. Крутизна потенциала подобрана из практических соображений. Этот парный потенциал (ниже – ЕАМ-1) показан на рис. 2. Особенность его – наличие справа от 1-го пика ПКФ (при 3.35 Å < r < 3.65 Å) участка функции с отрицательной кривизной, где вторая производная d²φ(r)/dr² < 0.

Потенциал погружения. Значения параметров эффективной электронной плотности ψ(r) были приняты равными p₁ = 4.4660, p₂ = 1.200, чтобы среднее значение $\left\langle \rho \right\rangle $ на атомах модели немного выше точки плавления было равно единице. Часть параметров потенциала погружения ЕАМ-1 была определена по температурной зависимости плотности вдоль бинодали, энергии и изотермическому модулю всестороннего сжатия металла (K_T). Модуль K_T сурьмы был рассчитан по значениям плотности [3] и скорости звука [5] с учетом МД-расчетов теплоемкостей C_V и C_p (соответственно 24.1 и 27.6 Дж/(моль К)); при 923 К он равен 20.6 ГПа. Согласно [4], при 933 К значение K_T = 21.6 ГПа. Параметры потенциала погружения ЕАМ-1, относящиеся к сжатым состояниям, определены по форме ударной адиабаты сурьмы (см. ниже). Все эти цифры приведены в табл. 3. Потенциал погружения сурьмы показан на рис. 3.

Расчетные свойства моделей сурьмы. МД-модели имели размер 2000 атомов в основном кубе. Релаксацию проводили алгоритмом Л. Верле с шагом по времени Δt = 0.01t₀, где t₀ = 1.123 × 10^{– 13} с. При расчетах учитывали тепловые вклады от электронов проводимости в энергию и давление по методу [29, 32], принимая число электронов проводимости на атом равным 5. Значения тепловой энергии электронов E_el приведены в табл. 4. Электронный вклад в давление p_el рассчитывали по формуле p_elV = (2/3)E_el [34]. Свойства моделей вдоль бинодали приведены в табл. 5. Получено хорошее согласие с опытом для плотности (d) и модуля K_T. При нагревании до 2000 К наблюдается постепенное занижение расчетной энергии модели сурьмы по сравнению с экспериментом [1, 2] (до 6.5 кДж/моль при 2023 К), частично обусловленное недостаточной точностью данных [1, 2] (теплоемкость принята равной 31.38 Дж/(моль К) при всех температурах до 1800 К).

В настоящей работе теплоемкости C_p и C_V рассчитывались методом МД с применением потенциала ЕАМ-1 на интервалах температур шириной 100 К, а модуль всестороннего сжатия – на интервалах объема шириной ∼0.5%. В табл. 6 приведены результаты расчетов теплоемкостей и скорости звука u_s (u_s = (K_TC_p/dC_V)^1/2) в моделях расплавов. В интервале 923–2023 К теплоемкости C_p и C_V меняются довольно слабо, а модуль K_T убывает в ~2 раза. Расчетная и фактическая скорости звука показаны на рис. 4. Несмотря на различную форму зависимости от температуры, расхождения с опытом по скорости звука не превышают ∼10%.

Плотность колебательных состояний (DOS) в модели сурьмы с потенциалом ЕАМ-1 при 1000 К монотонно убывает с ростом частоты ω (рис. 5), в отличие от расчетов [10], где получен промежуточный максимум DOS около энергии фонона 12 мэВ (при ωt₀ ≈ 2.05). Однако, DOS обращается в нуль и у нас, и в [10] при энергии фонона ~30 мэВ.

Коэффициенты самодиффузии. МД-расчеты показали, что потенциал ЕАМ-1 приводит к значению коэффициента самодиффузии D при 923 К (8.3 × 10^–5 см²/с), превышающему экспериментальное [7–9] в ∼1.5 раза (табл. 6). Завышение получается и при расчете через автокорреляционную функцию скоростей (9.80 × 10^–5 см²/с при 1000 К). Такое расхождение с опытом заслуживает внимания, так как обычно экспериментальные значения превосходят расчетные благодаря дополнительному вкладу от конвекции. С ростом температуры это расхождение уменьшается и при 1250 К полностью исчезает. Впрочем, экспериментальные значения D соседей сурьмы по Периодической системе также довольно высоки: D = = 7.72 × 10^–5 см²/с при 973 К для олова [36], 11.5 × × 10^–5 см²/с при 1073 К для висмута [37].

Используем довольно точное соотношение Стокса–Эйнштейна:

где η – динамическая вязкость и r_a – ионный радиус атома. Для жидких металлов величина r_a близка обычно к радиусу однозарядного иона [29]. В случае сурьмы при 923 К значение η = = 1.37 × 10^–3 Па с, и при D = 8.3 × 10^–5 см²/с мы получаем радиус r_a = 0.89 Å (в [6] было принято r_a ≈ 1.5 Å). Ионные радиусы для Sb³⁺ и Sb⁵⁺ равны соответственно 0.76 и 0.60 Å, а для иона Sb⁺ еще немного больше, так что значение D = 8.3 × × 10^‒5 см²/с оказывается вполне разумным. Отсюда следует, что данные [7, 8], возможно, требуют уточнения. Однако, завышение коэффициентов самодиффузии (вблизи от точки плавления) при использовании потенциала ЕАМ может быть обусловлено неучетом ковалентного вклада в межчастичное взаимодействие. Этот вклад занижает подвижность атомов, но его роль должна постепенно убывать с ростом температуры.

Структура моделей жидкой сурьмы. Согласно [6], побочный максимум ПКФ в районе 4.2 Å сохраняется при нагревании вдоль бинодали до ∼1100 К, далее до 1323 К виден в виде перегиба и исчезает при более высоких температурах. Как видно из рис. 1, расчетная и дифракционная ПКФ при 923 К практически совпадают. При 948 К наблюдается такое же совпадение двух ПКФ. Невязки R_g здесь очень малы (~0.035). При дальнейшем повышении температуры невязки остаются очень небольшими (0.039 при 1223 К и 0.047 при 1423 К), и согласие между дифракционными и модельными ПКФ остается очень хорошим вплоть до 1423 К (рис. 6). Наши ПКФ согласуются также с квантово-механическими расчетами [10] при всех температурах (до 1373 К).

Известно несколько способов расчета координационного числа (КЧ) в жидких металлах, однако для сурьмы, при наличии побочного максимума ПКФ, можно только применить вариант КЧ-2, когда КЧ равен удвоенной величине интеграла функции G(r) = 4πr²(N/V)g(r) на участке от r  = 0 до координаты (r_1g) первого максимума G(r). Найденные таким образом КЧ приведены в табл. 7. Значения КЧ невелики (3.56–3.94), что характерно для рыхлых систем типа жидких Si и Ge. Они очень слабо убывают с ростом температуры. Для сравнения, КЧ у металлов – соседей сурьмы по Периодической системе, рассчитанные тем же способом, значительно больше, чем у сурьмы: 8.85 (Sn) и 7.13 (Bi). Авторы [6] рассчитывали КЧ при 900–1400 К интегрированием до точки перегиба между первым пиком и побочным максимумом ПКФ и получили уменьшение КЧ от 6.6 до 6.1 с ростом температуры. Завышение по отношению к нашим значениям обусловлено асимметрией 1-го пика ПКФ. Метод расчета КЧ индивидуально для каждого атома [10] дает при 900–1300 К значения от 5.31 до 5.05.

Топологические характеристики модели сурьмы. Примем обозначения: r₁ – координата 1-го пика ПКФ, r_min – минимальное межчастичное расстояние в модели, r_mm – максимальное возможное значение r_min, равное 1.067(V/N)^1/3 [20, 38], g₁ = g(r₁) – высота 1-го пика ПКФ, r_1g – координата 1-го пика функции G(r) = 4πr²(N/V)g(r). Далее обозначим ρ₁ = r₁(N/V)^1/3, ρ_min = r_min(N/V)^1/3, ρ_mm = 1.067 и y = $\rho _{{\min }}^{{ - 1}} - r_{{{\text{mm}}}}^{{ - 1}}$. В [20, 38] показано, что для плотных некристаллических структур выполняется Первое топологическое правило, согласно которому, 1.067 < ρ₁ < 1.12. Для рыхлых структур (например, жидких и аморфных C, Si, Ge, As, Sb, S, Se, Te) ρ₁ < 1.02. Высота первого пика ПКФ определяется вторым правилом, согласно которому, выполняется соотношение:

Здесь g_1р – расчетное (“теоретическое”) значение высоты пика для плотных некристаллических структур. Для моделей аморфных структур с различными парными потенциалами отношение s = = g₁/g_1р = 1.035 ± 0.050, а для моделей жидкостей s = 1.001 ± 0.050. График уравнения (1) называется “главной последовательностью” [20, 38]. Точки для рыхлых структур типа непрерывной случайной сетки (S, Si, Te) располагаются выше главной последовательности.

В табл. 7 приведены характеристики структуры моделей жидкой сурьмы на бинодали. Величина ρ₁ плавно уменьшается при нагревании (от 0.967 до 0.931), отношение s также плавно убывает от 0.915 до 0.821. Только при переходе от 1175 до 1225 К наблюдается небольшой скачок s из-за смещения координаты r_min (от 2.443 к 2.398), которое связано скорее с методом обработки данных вблизи точки r_min. Таким образом, нет оснований делать вывод о структурных превращениях в модели сурьмы на бинодали. Однако низкие значения s < 1 (до 0.821) означают, что второе топологическое правило не выполняется для аномальных систем с побочным максимумом ПКФ. Переход части ближайших соседей атомов в отдельную “подгруппу субпика” приводит к понижению высоты 1-го пика ПКФ и к уменьшению отношения s.

Итак, межчастичный потенциал ЕАМ-1 хорошо описывает структуру жидкой сурьмы во всем исследованном интервале температур, и не требуется вводить дополнительных предположений об участии направленной связи (со своим специфическим потенциалом взаимодействия) и о повышении ригидности структуры с ростом температуры [6]. Аномальная форма ПКФ сурьмы (побочный максимум) объясняется отрицательной кривизной потенциала в интервале 3.35–3.65 Å (см. рис. 2), при которой межчастичные силы убывают при сближении частиц. Специфическое влияние отрицательной кривизны на свойства жидкостей обсуждается, например, в [39].

Эта особенность потенциала сурьмы не исчезает при нагревании. Она позволяет объяснить поведение аномального металла без привлечения механизмов, рассматриваемых для кристаллических систем с дальним порядком (Peierls distortion [6, 25]), и без участия ковалентной связи.

Построение моделей в условиях ударного сжатия. Характеристики ударной адиабаты сурьмы приведены в [40, 41]. Экспериментальная адиабата Гюгонио показана на рис. 7. По оси абсцисс отложена величина Z = V₀/V, где V – мольный объем сжатого металла, и исходный объем V₀ = = 18.178 см³/моль. Адиабата хорошо описывается выражением:

(критерий R² = 0.98733). По форме адиабаты были рассчитаны приведенные в табл. 3 оптимальные значения параметров потенциала погружения в ЕАМ-1, отвечающие за сжатые состояния. Эти параметры определили по методу [29, 32], исходя из условий, что давление (с учетом электронного вклада p_el) при данной степени сжатия должно равняться фактическому, а энергия модели (также с учетом электронного вклада E_el) должна равняться фактической, равной: Было учтено, что исходное давление p₀ ≈ 0. Кроме того, исходная энергия U₂₉₈ = –255.8 кДж/моль (табл. 5).

Свойства моделей сурьмы, построенных при условиях ударного сжатия с использованием потенциала ЕАМ-1, приведены в табл. 8 и на рис. 7, 8. В целом получено очень хорошее согласие с опытом по энергии (колонки 8 и 9) и давлению (колонки 2 и 11). Температура и коэффициент самодиффузии плавно возрастают вдоль ударной адиабаты (см. рис. 8).

Структура под давлением. В работе [14] были исследованы структурные свойства жидкой сурьмы под давлением до 9.4 ГПа при температурах на 50 К выше линии плавления методом дифракции синхротронного излучения. Они не вполне согласуются с данными [6] при обычном давлении, где координата 1-го пика вблизи T_m близка к 3.04 Å; в [14] эта координата близка к 2.98 Å. Соответственно координата 2-го пика равна 6.35 Å в [6] и 6.40 в [14]. Можно сравнить наши расчеты с данными [14]. Мы построили методом МД с потенциалом ЕАМ-1 модели сурьмы при давлениях 3.0, 4.1, 6.5 и 9.4 ГПа и при температурах соответственно 925, 855, 930 и 1030 К (на 50 К выше линии плавления, согласно [42]). Электронные вклады в давление при этих температурах невелики (до 0.032 ГПа при 1030 К). При давлении 9.4 ГПа плотность модели составила 8.059 г/см³. Координаты пиков ПКФ этой модели равны r₁ = = 2.995 Å и r₂ = 6.01 Å, что немного больше данных [14]: r₁ = 2.94 ± 0.01 и r₂ = 5.95 ± 0.03 Å. Однако имеется значительное расхождение по форме графиков дифракционных и модельных ПКФ по высоте 1-го пика (см. рис. 9). Высота 1-го пика ПКФ модели монотонно растет с давлением от 2.38 до 2.56, а высота пиков в [14] практически постоянна (∼2.0). Невязки по ПКФ здесь довольно велики: от R_g = 0.118 при давлении 3.0 ГПа до 0.166 при давлении 9.4 ГПа.

Координационные числа модели рассчитывали по методу КЧ-2 (см. выше). При указанных выше давлениях КЧ = 4.76, 5.14, 5.37 и 5.72. Авторы [14] получили КЧ = 4.7 при 0.7 ГПа и 5–7 при более высоких давлениях. Эти данные в общем согласуются, несмотря на расхождения в форме ПКФ. Однако в целом результаты моделирования значительно отличаются от данных дифракционного исследования [14]. Для выяснения причин этих расхождений была дополнительно построена серия моделей сурьмы на изотерме 930 К при давлениях от 0 до 8 ГПа (потенциал ЕАМ-1). Некоторые свойства этих моделей приведены в табл. 9. Значения КЧ (метод КЧ-2) постепенно увеличиваются от 3.92 до 5.72 с ростом плотности, а приведенная координата первого пика ρ₁ возрастает при этом от 0.967 до 1.023 (несмотря на убывание r₁, в отличие от данных [14], где r₁ увеличивается от 2.98 до 3.01 Ǻ с ростом давления до 5 ГПа). Судя по этим величинам, структура жидкой сурьмы рыхлая при всех давлениях до 9.4 ГПа, но степень рыхлости значительно меньше, чем на бинодали. Кроме того, высота 1-го пика ПКФ модели увеличивается с ростом давления (в [14] от давления не зависит), а приведенная высота 1-го пика s всюду близка к единице. Это указывает на хорошую выполнимость второго топологического правила у моделей и на отсутствие признаков ковалентной связи, что не согласуется с концепцией [14]. Мы не видим также никаких особенностей вблизи 3 и 5 ГПа, которые указывали бы, согласно [14], на значительную роль ковалентной связи в этих условиях.

Расчеты по дифракционным данным [11, 12]. На дифракционной ПКФ сурьмы при 933 К, полученной в [11, 12], также наблюдается побочный максимум вблизи 4.8 Ả. Алгоритм Шоммерса приводит к модельной ПКФ с невязкой R_g = = 0.049. Назовем полученный потенциал ЕАМ-2. Кривизна этого потенциала отрицательна в интервалах 2.80–2.85, 3.70–3.90, 4.45–4.55 и 5.05–6.25 Ǻ. При релаксации модели сурьмы этот вариант дает коэффициент самодиффузии при 933 К, равный 4.78 × 10^–5 см²/с, что близко к опытным данным (4.45–5.65) × 10^{– 5} см²/с при 933 К [7, 8]. И здесь причиной аномалии является отрицательная кривизна потенциала. Однако при переходе от ЕАМ-1 к ЕАМ-2 “выигрыш” в коэффициентах самодиффузии сопровождается проигрышем в отношении координаты 1-го пика, значительно отличающейся от всех остальных структурных исследований (см. табл. 1). Поэтому вариант ЕАМ-2 мы считаем пока менее вероятным.

Таким образом, на основе модели погруженного атома (ЕАМ) можно удовлетворительно описать большинство свойств жидкой сурьмы с рыхлой структурой как вдоль бинодали, так и в условиях ударного сжатия (плотность, энергия, модуль всестороннего сжатия, скорость звука, ПКФ, координационные числа, ударная адиабата). Предложенный выше потенциал ЕАМ-1 позволяет рассчитать свойства жидкой сурьмы в широком интервале температур и давлений (эти данные не представлены в статье только из-за недостатка места). Оказалось возможным использовать центрально-симметричный потенциал ЕАМ-1 без включения направленной связи и соображений, заимствованных из зонной теории (Peierls distortion). Причиной аномалий структуры сурьмы является наличие интервалов расстояний, на которых межчастичный потенциал имеет отрицательную кривизну. Однако потенциал ЕАМ-1 не приводит к согласию с дифракционными данными [14] в области давлений до 10 ГПа и с измерениями коэффициента самодиффузии в [7]. Причины этих расхождений пока неясны, и необходимы новые экспериментальные данные.

Автор благодарит проф. Г. Макова (Израиль) за предоставление таблиц ПКФ жидкой сурьмы при 923 и 948 К (ПКФ при остальных температурах были получены оцифровкой графиков из [6]).