Журнал физической химии, 2021, T. 95, № 2, стр. 294-296

Распад метастабильного состояния начинается с образования и последующего роста зародыша новой фазы. Работа образования жизнеспособного (критического) зародыша определяется высотой активационного барьера, отделяющего метастабильную фазу от стабильной, и может быть рассчитана в рамках теорий капиллярности Гиббса [1] и Ван-дер-Ваальса [2], метода функционала плотности [3], при молекулярно-динамическом моделировании [4–10]. В последнем случае используются методы зонтичной выборки [4, 5], выбора путей перехода [6, 7], внедрения зародыша [8], метадинамики [9, 10] и другие.

В случае кристаллического зародыша равновесная форма его огранки определяется правилом Кюри–Вульфа [11]. Согласно теории капиллярности Гиббса, работа образования критического зародыша [12]:

где $\gamma _{e}^{{}}$ – эффективная поверхностная свободная энергия, ${{p}_{*}}$ – давление в критическом зародыше, ${{p}_{l}}$ – давление в метастабильной жидкости, ${{R}_{*}}$ – радиус критического зародыша. Уравнение (1) задает сферическую аппроксимацию кристаллического зародыша.

Работа образования критического зародыша может быть также определена из основного уравнения классической теории нуклеации (КТН):

где J – частота зародышеобразования, ${{\rho }_{l}}$ – числовая плотность жидкости, B – кинетический коэффициент, ${{k}_{{\text{B}}}}$ – постоянная Больцмана, T – температура.

В данной работе методом молекулярной динамики (МД) исследуется кристаллизация переохлажденной леннард-джонсовской (ЛД) жидкости. Для определения параметров кристаллизационного процесса привлекается метадинамика, которая ускоряет переход системы из метастабильного состояния в стабильное и позволяет определить рельеф термодинамического потенциала в пространстве обобщенных координат задачи.

Исследуемая система содержала N = 5324 ЛД-частиц, которые размещались в кубической ячейке с периодическими граничными условиями. Параметры потенциала σ и ε, масса частицы m и постоянная ${{k}_{{\text{B}}}}$ использовались для приведения рассчитываемых величин к безразмерному виду. Далее все величины представляются в приведенном виде. Радиус обрезания межчастичных взаимодействий принимался равным r_c = 6.78. Расчеты проводились в NpT-ансамбле по изобарам ${{p}_{l}}$ = = –0.995, 0 и 6.912. Для выбранных давлений температура равновесного фазового перехода ${{T}_{m}}$ = = 0.6, 0.692 и 1.2. На изобарах ${{p}_{l}}$ = –0.995 и 0 кристаллизация исследовалась при температуре ниже температуры конечной точки линии плавления – точки встречи ее метастабильного продолжения со спинодалью перегретой жидкости ${{T}_{K}}$ (${{T}_{K}} = 0.529$) [13]. Использовался пакет параллельного МД-моделирования LAMMPS [14]. Шаг интегрирования уравнений движения частиц по времени составлял Δt = 0.002318. Для реализации метадинамического процесса применялась библиотека открытого доступа PLUMED версии 2.6.0 [15].

Метадинамика увеличивает вероятность наблюдения состояний с высокими значениями термодинамического потенциала. Это достигается добавлением к гамильтониану системы зависящего от времени потенциала смещения $V[\vec {s}({{\vec {r}}^{N}},t)]$ – функции коллективных переменных, задаваемых векторной функцией $\vec {s}({{\vec {r}}^{N}})$ положений N частиц системы [9]. Введение потенциала смещения приводит к снижению вероятности повторного исследования уже посещенных областей коллективных переменных. В пределе больших времен потенциал смещения “проталкивает” систему над активационным барьером. Термодинамический потенциал при этом $\Phi (\vec {s}) = - \mathop {\lim V}\limits_{t \to \infty } [\vec {s}({{\vec {r}}^{N}},t)]$.

В качестве коллективных переменных использованы средний параметр порядка Стейнхардта Q₆ [16], рассчитываемый по 350 частицам вокруг случайно выбранной частицы, и потенциальная энергия системы U, которая вычислялась в процессе МД-моделирования. Потенциал $V[\vec {s}({{\vec {r}}^{N}},t)]$ задавался как ряд гауссовских функций, вносимых с шагом τ по времени. Временной шаг добавления гауссианов 0.023. Ширина гауссианов 0.01 для параметра Стейнхардта и 40 для потенциальной энергии, высота гауссианов $1{\text{ }}{{k}_{{\text{B}}}}T$.

При заданных значениях ${{p}_{l}}_{{}}$ и T  регистрировалось до 40 событий кристаллизации жидкости, по которым вычислялось среднее значение высоты активационного барьера. Для всех исследованных состояний определены поверхности термодинамического потенциала в пространстве коллективных переменных Q₆, U. Переход системы из метастабильной потенциальной ямы в более глубокую потенциальную яму твердого тела связан вначале с совершенствованием порядка в ограниченном объеме (увеличение Q₆ при слабом изменении U), а затем, после прохождения переходного состояния, c кристаллизацией всей системы (резкое снижение U). Высота активационного барьера кристаллизации $\Delta \Phi $, равная работе образования критического зародыша ${{W}_{*}}$, определяется как разность между максимальным значением термодинамического потенциала $\Phi $ в переходном состоянии и его минимальным значением, отвечающим дну потенциальной ямы метастабильного состояния.

Зависимость приведенной работы образования критического зародыша ${{W}_{*}}{\text{/}}T$ от температуры представлена на рис. 1. Там же полученные данные сопоставляются с расчетами ${{W}_{*}}{\text{/}}T$ из КТН (уравнение (2)) по МД-данным о частоте нуклеации [17] и в рамках теории капиллярности Гиббса (уравнение (1)) без учета зависимости поверхностной свободной энергии критического зародыша ${{\gamma }_{*}}$ от его размера ${{R}_{*}}$, т.е. в приближении ${{\gamma }_{*}} = {{\gamma }_{\infty }}$. Значения поверхностной свободной энергии ЛД-системы на плоской межфазной границе ${{\gamma }_{\infty }}$ при температуре выше температуры конечной точки линии плавления получены в работе [18]. Ниже этой температуры сосуществование жидкости и кристалла на плоской межфазной границе становится невозможным [13, 19]. Давление внутри критического зародыша определялось из условия равенства химических потенциалов переохлажденной жидкости и кристаллического зародыша. Термическое уравнение состояния ЛД-системы, а также алгоритм расчета химических потенциалов представлены в работах [13, 19].

Как следует из рис. 1, полученные с помощью метадинамики значения ${{W}_{*}}{\text{/}}T$ хорошо согласуются с рассчитанными из КТН по молекулярно-динамическим данным о частоте кристаллизации. В то же время имеет место систематическое рассогласование результатов метадинамики и теории капиллярности Гиббса в приближении ${{\gamma }_{*}} = {{\gamma }_{\infty }}$, т.е. без учета размерной зависимости поверхностной свободной энергии критических зародышей. Величина рассогласования растет по мере повышения давления и температуры системы.

Метадинамические данные о ${{W}_{*}}$ и ${{p}_{*}}$ позволяют рассчитать радиус ${{R}_{*}}$ и эффективную поверхностную свободную энергию критического кристаллического зародыша ${{\gamma }_{e}}$ (уравнение (1)). Зависимость ${{\gamma }_{e}}$ от кривизны межфазной границы, при условии постоянства давления в жидкой фазе, показана на рис. 2. Расхождения между ${{\gamma }_{e}}$ и ${{\gamma }_{\infty }}$ увеличиваются с ростом давления. При ${{p}_{l}} = 6.912$, когда ${{R}_{*}} \approx 3.0$, величина ${{\gamma }_{e}}$ на 20% превышает ${{\gamma }_{\infty }}$.

Таким образом, рост с повышением давления расхождения данных по МД-моделированию нуклеации в переохлажденной жидкости и теории капиллярности Гиббса связан с “размерным эффектом” − зависимостью поверхностной свободной энергии кристаллических зародышей от их размера. Метадинамика может быть эффективным методом исследования размерного эффекта при малых переохлаждениях, где невозможен прямой молекулярно-динамический расчет частоты нуклеации.

Работа выполнена при финансовой поддержке проекта РФФИ-Урал № 20-48-660027 р_а. Расчеты проведены с использованием суперкомпьютера “Уран” ИММ УрО РАН.