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

МЕТОДИКА АНАЛИЗА

Комбинированный метод ФМ + КХП

Рассматривается система, состоящая из N_fr-узлов фрагмента и одного виртуального распределенного узла, заполнение которого описывается с помощью КХП. Каждый узел фрагмента может быть свободен или занят частицей. Состояние заполнений всех узлов фрагмента рассматривается как одна квазичастица, имеющая конфигурацию k. Это состояние задается функцией $\sigma (k,f)$ равной нулю или единице, если узел f в данной конфигурации k свободен или занят, где $1 \leqslant f \leqslant {{N}_{{fr}}}$ [18]. Всего имеется разных $B = {{2}^{{{{N}_{{fr}}}}}}$, 1 ≤ k ≤ B, конфигураций. Вероятность того, чтобы найти k-ю конфигурацию частиц на фрагменте принимает форму [12–14]

$\begin{gathered} {{\theta }_{k}} = \exp ( - \beta {{E}_{k}}){\text{/}}Q, \\ Q = \sum\limits_{k = 1}^B {\exp ( - \beta {{E}_{k}})} , \\ \end{gathered} $

здесь сумма по  f, перечисляя все узлы фрагмента, определяет конфигурацию k, Q_f – энергия связи частицы с решеткой в узле  f, E_k – полная энергия молекул в данной конфигурации, функция $\sigma (k,f)$ равна нулю или единице, если узел  f, 1 ≤ ≤ f ≤ ${{N}_{{fr}}}$, в данной конфигурации k свободен или занят;  f и h нумеруют узлы фрагмента, второе слагаемое есть энергия латеральных взаимодействий частиц в конфигурации k (отсутствие частиц в соседних узлах f и h дает нулевой вклад); штрих в сумме по h исключает узел f; символ g обозначает распределенный узел, описываемый в КХП, ${{z}_{{fg}}}$ – число связей между узлом f и узлом g, зацепление заселенности узлов ФМ и КХП реализуется через условные вероятности для парных функций t_fg; Q – статистическая сумма системы.

Концепция квазичастиц подразумевает, что смесь квазичастиц идеальна, и все термодинамические и кинетические функции вычислены уравнениями для идеальных смесей [18]. Концентрация квазичастицы определяется как

${{\theta }_{k}} = \sum\limits_{n = 1}^B {\exp \{ \beta [{{E}_{n}} - {{E}_{k}}]\} } ,\quad \sum\limits_{k = 1}^B {{{\theta }_{k}}} = 1.$

Все дальнейшие операции усреднения по конфигурациям выражаются через вероятности θ_k.

Локальная плотность узла θ_f , 1 ≤ f ≤ ${{N}_{{fr}}}$ и средняя плотность фрагмента θ_fr запишутся как ${{\theta }_{f}} = \sum\nolimits_{l = 1}^B {\sigma _{{f,l}}^{{}}{{\theta }_{l}}} $ и ${{\theta }_{{fr}}} = \sum\nolimits_k {{{n}_{k}}{{\theta }_{k}}{\text{/}}{{N}_{{fr}}}} $, где через n_k обозначено число занятых узлов на фрагменте. Формулу для средней плотности фрагмента θ_fr можно также записать как ${{\theta }_{{fr}}} = \sum\nolimits_{f = 1}^{{{N}_{{fr}}}} {{{\theta }_{f}}{\text{/}}{{N}_{{fr}}}} $, тогда средняя плотность системы выражается по формуле для обычных неоднородных систем [11, 12] в виде

${{\theta }_{{sys}}} = \sum\limits_{f = 1}^{{{N}_{{fr}}}} {{{\theta }_{f}}{\text{/}}N + {{\theta }_{g}}{\text{/}}N} ,\quad N = {{N}_{{fr}}} + 1.$

Эта запись отражает заполнение каждого узла системы. Выражение локальной изотермы для узла g имеют вид:

${{a}_{0}}P = \left( {\frac{{{{\theta }_{g}}}}{{1 - {{\theta }_{g}}}}} \right){{\Lambda }_{g}},\quad {{\Lambda }_{g}} = \prod\limits_{f \in {{z}_{g}}} {{{S}_{{gf}}}} ,$

(1)

$\begin{gathered} {{S}_{{gf}}} = 1 + {{t}_{{gf}}}{{x}_{{gf}}},\quad {{x}_{{gf}}} = \exp ( - \beta {{\varepsilon }_{{gf}}}) - 1, \\ {{t}_{{gf}}} \equiv t_{{gf}}^{{AA}} = \theta _{{gf}}^{{AA}}{\text{/}}\theta _{g}^{A} = 2\theta _{f}^{A}{\text{/}}[{{\delta }_{{gf}}} + {{b}_{{gf}}}], \\ \end{gathered} $

$\begin{gathered} {{\delta }_{{gf}}} = 1 + {{x}_{{gf}}}(1 - {{\theta }_{f}} - {{\theta }_{g}}), \\ {{b}_{{gf}}} = {{[{{\delta }_{{gf}}}^{2} + 4{{x}_{{gf}}}{{\theta }_{f}}{{\theta }_{g}}]}^{{1/2}}}, \\ \end{gathered} $

где θ_g ≡ $\theta _{g}^{A}$ – вероятность заполнения узла g частицей А, $\theta _{g}^{V} = 1 - {{\theta }_{g}}$ – локальная доля вакансий; $\theta _{{gf}}^{{in}}$ – вероятность нахождения частицы i на узле g и частицы n на узле f; ${{t}_{{gf}}}$ – условная вероятность нахождения частицы А рядом с другой частицей А. Здесь символ соседнего узла f пробегает всех соседей z_g узла g. В формуле (8) имеет место соотношение ${{t}_{{gf}}} = \theta _{{gf}}^{{AA}}{\text{/}}\theta _{g}^{A} = \theta _{f}^{A}{{t}_{{fg}}}{\text{/}}\theta _{g}^{A},$ следующее из естественной симметрии $\theta _{{gf}}^{{AA}} = \theta _{{fg}}^{{AA}}$, что определяет связь между заполнениями узлов  f  фрагмента и КХП g.

Рассматриваемая делокализация узла КХП в случае однородного фрагмента означает, что все локальные плотности фрагмента одинаковые, как и в случае отсутствия данного узла КХП. Поэтому введение такого узла можно считать как введение калибровочной функции, которая позволяет повысить точность описания распределения частиц в широком диапазоне температур и плотностей [10]. Трактовка распределенного узла КХП в качестве калибровочной функции для уравнений совместного применения ФМ + КХП позволяет найти ее параметры по известному точному значению критической температуры. Каждый узел f фрагмента имеет дополнительный вклад в его энергию от узла КХП g. Это вводится через одночастичный вклад в энергию узла – добавляется величина равная ${{\varepsilon }_{{fg}}}{{z}_{{fg}}}{{t}_{{fg}}}$, где энергетический параметр ε имеет свое обычное значение, а число “соседей” z_fg можно определить в виде z_fg = K/N_fr , где K – некоторый числовой параметр. Если параметр K имеет смысл числа соседей узла КХП (K = z), то z_fg = z/N_fr . В общем случае число K можно ассоциировать с эффективным числом соседей отличным от рассматриваемой решетки. Очевидно, если K = 2, то мы имеем аналог одномерной системы, в которой отсутствуют фазовые переходы. При K = 0 вклады КХП-узлов отсутствуют и мы получаем переход к выражениям для исходного ФМ. С учетом этого, выражение для локальной изотермы единственного виртуального узла g выписывается в прежнем виде уравнения (1), в котором функция неидеальности запишется как ${{\Lambda }_{g}} = \prod\nolimits_{f \in {{N}_{{fr}}}} {{{{({{S}_{{gf}}})}}^{{{{z}_{{fg}}}}}}} = {{({{S}_{{gf}}})}^{K}},$ поэтому выражение (1) упрощается. В результате параметрами системы уравнений становятся параметры уравнения КХП (1) – указанный выше параметр К (для числа соседей) и величина энергии связи частицы на узле КХП с решеткой Q_g.

Алгоритм и условия расчета

При проведении расчетов функция $\sigma \left( {k,f} \right)$ задавалась в виде массива нулей и единиц по формуле $\sigma \left( {k,f} \right)$ = остаток от деления $\operatorname{int} \{ (k - 1){\text{/}}{{2}^{{f - 1}}}\} $ на 2, где f номер узла фрагмента [19]. При этом k‑я строчка массива определяет состояние занятости узла в конфигурации k. Например, для N_fr = 4 пятая строка массива $\sigma \left( {k,f} \right)$ имеет вид (0010), что означает наличие в конфигурации k = 5 занятого частицей узла   f = 3, а остальные узлы – свободны.

Для ускорения расчетов вместо указанной формулы для $\sigma \left( {k,f} \right)$ можно пользоваться двоичным представлением десятичных чисел (k – 1), 1 ≤ k ≤ B, в котором каждая цифра 0 или 1 с порядковым номером  f, отсчитываемым справа налево, задает состояние $\sigma \left( {k,f} \right)$. Если в первом случае нам для каждого узла f нужно использовать формулу с выделением остатка, то во втором случае одна операция переводит десятичное число (k – 1) в бинарную строчку, которая содержит ту же информацию, что формула (это программа сделает быстрее). Например, для N_fr = 4 и конфигурации k = 5 двоичное представление десятичного числа (k – 1) = 4 имеет вид (0100), где цифры нумеруем справа налево. Тогда мы также получаем, что в конфигурации k = 5 частицей занят узел f = = 3, а остальные свободны. Преимущество этого описания заключается в том, что можно без перебора узлов  f  осуществлять операции над значениями $\sigma \left( {k,f} \right)$ сразу для всех узлов f для заданной конфигурации k посредством операций над двоичными числами, являющихся представлением десятичных (k – 1). Метод позволил без больших затрат времени рассчитать фрагменты до N_fr = 30 на обычном персональном компьюторе.

Использованные в работе энергии отсчитаны от связанного состояния атомов, т.е. притяжению отвечает положительный знак, а отталкиванию отрицательный. В представленных ниже результатах потенциал решетки в узлах систем Q_f  отсчитывался либо от 0, либо от βQ_f/ε = 22 (что формально отвечает изотермам адсорбции). В последнем случаем изотермы имеют параллельный сдвиг по давлению к меньшим величинам по сравнению с Q_f = 0, но фазовые диаграммы не зависят от этого сдвига по давлениям.

РЕЗУЛЬТАТЫ ЧИСЛЕННОГО АНАЛИЗА

Надкритические температуры. ФМ не позволяет получить расслаивание фаз, поэтому простейшая ситуация для анализа влияния непрямых корреляций соответствует надкритической области температур. На рис. 1 приведены кривые изотерм, а также корреляторы порядка n = 2, 3 и 5 полученные в разных приближениях ФМ и КХП для квадратной решетки z = 4, d = 2 при βQ_f/ε = 22.

Рис. 1.

Изотермы плотности (а) и корреляторов второго (б), третьего (в) и пятого (г) порядков. Обозначения см. текст.

Изотермы на рис. 1а представлены пунктирными кривыми по расчетам в КХП (1, 6) и сплошными кривыми ФМ для фрагментов 2 × 2 (2 и 7), 3 × 3 (3 и 8), 4 × 4 (4 и 9) и 5 × 5 (5 и 10) в широкой области флюида от разреженного пара до плотного флюида для температуры τ = 1.3 (1–5), τ = = T/T_c(КХП), T_c(КХП) – критическая температура в КХП. (Для сопоставления также приведены изотермы в подкритической области при τ = 0.7 (6–10)).

Из рис. 1а видно, что увеличение фрагмента в размере меняет угол наклона γ изотермы в точке θ = 0.5 при неизменном давлении в ней. Для количественного сопоставления изотерм в разных способах расчета воспользуемся величиной тангенса угла наклона изотермы $\operatorname{tg} \gamma $ ≡ $\operatorname{tg} {{(\gamma )}_{{\theta = 0.5}}}$ при плотности θ = 0.5, отчет угла γ идет от наклона к оси плотности. Для кривых на рис. 1а имеем значения $\operatorname{tg} \gamma $ показанные в табл. 1.

С ростом размера фрагмента при температуре τ = 1.3 угол наклона изотермы уменьшается и создается впечатление, что он приближается к углу наклона для КХП. Однако эта тенденция обусловлена двумя разными факторами. Увеличение размера фрагмента увеличивает только учет эффектов корреляции, что приводит к уменьшению угла, а для КХП данная область температур уже достаточно близка к критической температуре в данном приближении, что приводит к существенному уменьшение угла вблизи Т_c (КХП). Влияние температуры при расчете в ФМ для фрагментов 3 × 3 и 4 × 4 и в КХП на угол наклона показан в табл. 2.

При высоких температурах три кривые дают достаточно близкие значения углов наклона изотермы в КХП и ФМ 3 × 3 и 4 × 4. Причем $\operatorname{tg} \gamma $ ФМ 3 × 3 > $\operatorname{tg} \gamma $ ФМ 4 × 4 > $\operatorname{tg} \gamma $ КХП. Данная тенденция сохраняется с уменьшением температуры, пока изотерма не входит в область влияния критической точки. Расхождение между $\operatorname{tg} \gamma $ (ФМ 4 × × 4) > $\operatorname{tg} \gamma $ КХП резко увеличивается в области от τ = 1.5 в ходе резкого уменьшения $\operatorname{tg} \gamma $ КХП.

В то же время следует учитывать, что и ФМ является вблизи критической области приближенным методом, в котором могут проявляться дальние корреляции. Это демонстрируют данные для τ = 0.7 в табл. 1. В подкритической области углы наклонов достаточно заметно отличаются от вертикального положения секущей и, тем более, от отрицательного угла наклона изотремы. Здесь даже расчет для фрагмента 5 × 6 приводит к углу $\operatorname{tg} \gamma $ (ФМ 5 × 6) = 0.16.

Более детально отличие эффектов непрямых корреляций от прямых корреляций в КХП представляют собой отношения корреляторов разной размерности. Корреляторы порядка n определяются как средние величины по фрагменту ${{\theta }_{{(n)}}} = \sum\nolimits_{k = 1}^B {N_{{(n)}}^{A}(k){{\theta }_{k}}} {\text{/}}{{N}_{{(n)}}},$ где $N_{{\left( n \right)}}^{A}(k)$ – число групп из n соседних занятых узлов фрагмента для его конфигурации k, ${{N}_{{(n)}}}$ – число групп из n числа соседних узлов во фрагменте. Так как фрагмент однороден, то достаточно рассчитывать корреляторы, построенные на одном центральном узле и на его (n–1) узлах окружения. Когда имеются разные способы расположений соседних узлов вокруг центра (неприводимые конфигурации [11]), то здесь по ним проводится дополнительное усреднение. Так для тройного коррелятора ${{\theta }_{{(3)}}} = \theta _{{fG(f)}}^{{AAA}} = \sum\nolimits_{k = 1}^B {N_{{fG(f)}}^{{AAA}}\left( k \right){{\theta }_{k}}} {\text{/}}6,$ где функция рассчитывается на конкретном центральном узле f   и его полном окружении G(f). То есть рассматривается полное окружение центрального узла с двумя соседними узлами, чтобы учесть разное расположение трех соседних узлов системы. Всего вариантов комбинации трех узлов 6, поэтому в знаменателе 6.

Роль размера фрагмента и выявление роли непрямых корреляций в изотермах демонстрируют величины ξ_n, представляемые в нормированном виде через отношения корреляторов ξ_n = = θ_(n)(ФМ)/θ_(n)(КХП) в сопоставляемых приближениях ФМ и КХП.

На рис. 1б,в,г представлены для ФМ изотермы корреляторов парного (рис. 1б), третьего (рис. 1в), и пятого (рис. 1г) порядков, нормированных на соответствующие значения корреляторов в КХП, при τ = 2.5 (1 и 5), 2.1 (2 и 6), 1.7 (3 и 7) и 1.3 (4 и 8) для фрагментов 3 × 3 (сплошные линии 1–4) и 4 × 4 (пунктирные линии 5–8). При высоких температурах τ = 2.5 (1 и 5) и 2.1 (2 и 6) отношение корреляторов ξ_n монотонно возрастает с ростом плотности, стремясь к единице при полном заполнении системы, когда результаты по КХП и ФМ совпадают. С убыванием температуры отношения ξ_n растут и при низких температурах τ = 1.7 (3 и 7) и 1.3 (4 и 8) корреляторы изменяются с образованием максимума, сохраняя стремление к единице при больших плотностях, но уже стремление сверху. Чем выше порядок корреляторов, тем выше максимум значения их отношения: на рис. 1г максимум самый большой. Это показывает, как с ростом порядка корреляторов растет расхождение между ФМ и КХП. Данное различие еще больше возрастает при стремлении к критической изотерме сверху (например, при τ = 1.1 максимум на поле 1г будет порядка 32 вместо 8!).

Влияние параметров калибровочных функций

Совместное использование ФМ + КХП позволяет получить петлю на изотерме, которая отвечает кривым типа ван-дер-Ваальса, и из анализа петли можно определить сосуществующие плотности разреженной и плотной фаз. Зная сосуществующие плотности при каждой температуре, построим кривую расслаивания в рамках обсуждаемой модели. Уравнения модели содержат два молекулярных параметра, относящиеся к характеристкам распределенного узла КХП, которые фактически являются параметрами калибровочной функции: величина K (для числа соседей) и величина энергии связи частицы на узле КХП с решеткой Q_g (1). На рис. 2 приведены примеры влияния этих параметров на изотермы (поля 2а и 2в) и кривые расслаивания (поля 2б и 2г) для фрагмента 3 × 3. На всех полях приведены аналогичные значения в КХП (пунктирная линия) и аналогичные изотермы в ФМ (точечная линия).

Рис. 2.

Изотермы (а, в) и фазовые диаграммы (б, г) в ФМ + КХП.

На рис. 2а построены изотермы плотности фрагмента (1–3) (сплошные линии) и плотности виртуального узла (4–6) (пунктирные линии) для ФМ + КХП на фрагменте 3 × 3, полученные при K = 4 (1,4), 3.5 (2, 5) и 3 (3, 6), а также изотерма с секущей в КХП (пунктирная кривая 7) и изотерма ФМ (точечная кривая 8) при βQ_f/ε = 22. С уменьшением K уменьшается к.ч. виртуального узла, поэтому уменьшаются его заполнения.

На поле 2б показаны соответствующие фазовые диаграммы для ФМ + КХП на фрагменте 3 × × 3 при тех же параметрах K для фргамента 3 × 3 и в КХП (пунктир). Локальные плотности узлов фрагмента и КХП отличаются между собой за счет различий в учете непрямых корреляций. В целом, хотя отличия между локальными плотностями не очень велики, кривые рассливания отличаются по Т_с значительно. Отметим, приближенный характер получаемых кривых расслаивания (кривая 3) может иметь темпетатуру выше точного решения Онсагера, которой соответствует величина $(\beta \varepsilon )_{{\text{c}}}^{{ - 1}}$ = 0.567. Поэтому кривая (3), имеющая температуру ниже данной величины, отвечает нефизическому значению.

На рис. 2в показаны изотермы плотности фрагмента (1–4) при Q_f   = 0 и плотности виртуального узла (5–8) для ФМ + КХП с потенциальной энергией для виртуального узла Q_g = 0 (1, 5), 1/3ε (2, 6), –1/3ε (3, 7) и –ε (4, 8), а также КХП (пунктирная кривая с секущей 9) и ФМ (точечная кривая 10), при K = 4.

Локальные плотности узлов фрагмента и виртуального узла отличаются между собой за счет различий в учете непрямых корреляций. В варианте Q_g = 1/3ε (2, 6) локальные изотермы для фрагмента и виртуального узла ближе всего между собой. Тогда как для вариантов Q_g = 0 (1, 5) и 2/3ε (3, 7) эти отличия более заметны. Здесь следует учитывать, что вклад от виртуального узла быстро уменьшается с ростом размера фрагмента.

На рис. 2г показаны фазовые диаграммы для ФМ + КХП с потенциальной энергией для виртуального узла Q_g = 0 (1), 1/3ε (2), –1/3ε (3) и –ε (4), а также КХП (пунктирная кривая 5). Увеличение притяжения Q_g > 0 смещает кривую расслаивания влево и увеличивает ее асимметрию за счет уменьшения величины критического заполнения ${{\theta }_{c}}$. Отталкивание между решеткой и частицей Q_g < 0 увеличивает ${{\theta }_{c}}$.

Квазиточные кривые расслаивания

На рис. 3а–в показаны кривые расслаивания для трех фрагментов 2 × 2 (рис. 3а), 3 × 3 (рис. 3б) и 4 × 4 (рис. 3в), сопоставляемые с расслаиванием в КХП (пунктирные линии).

Рис. 3.

Фазовые диаграммы (а–в) и изотермы (г) в ФМ + КХП с калибровкой.

На рис. 3а–в кривые расслаивания 1 отвечают K = 4, Q_g = 0. С увеличением размера фрагмента величина Т_c уменьшается, приближаясь к точному значению: фрагмент 2 × 2 (βε)_c = 1.536, θ_c = = 0.396, фрагмент 3 × 3 (βε)_c = 1.642, θ_c = 0.368, фрагмент 4 × 4 (βε)_c = 1.712, θ_c = 0.355.

Кривые 2 отвечают значениям K = 3.07, θ_c = = 0.358 (2 × 2), K = 3.39, θ_c = 0.342 (3 × 3) и K = = 3.69, θ_c = 0.343 (4 × 4) калиброванным под точное значение (βε)_c = 1.764 при Q_g = 0. С ростом размера фрагмента величина параметра К стремится к своему значению K = z = 4.

Кривые 3 отвечают калиброванным значениям K и Q_g под точное значение (βε)_c = 1.764 и под значение критической плотности θ_c = 0.500 для всех кривых (чтобы сохранить тенденцию к сохранению симметрии кривой расслаивания, которая нарушается в общем случае в ФМ + КХП). Это дает для фрагмента 2 × 2: K = 3.79, Q_g = –1.38, для фрагмента 3 × 3: K = 3.61, Q_g = –1.23, и для фрагмента 4 × 4: K = 3.72, Q_g = –1.19. С ростом размера фрагмента величина K слабо меняется, а величина Q_g  уменьшается по модулю.

На рис. 3г представлены изотермы для фрагментов 2 × 2 (1, 4), 3 × 3 (2, 5) и 4 × 4 (3, 6) с калибровкой в ФМ + КХП (1–3), без калибровки в ФМ (4–6) и в КХП (пунктирная кривая с расслаиванием 7) при Q_f = 0. Критические изотермы ФМ + + КХП (1–3) построены при калибровочных параметрах K и Q_g, полученных для кривых 3 на рис. 3а–в, обеспечивающих точное значение (βε)_c = = 1.764 и критическую плотность θ_c = 0.500. С ростом фрагмента для ФМ + КХП получен рост давления: кривая 1 лежит левее всех в области малых давлений, следующие кривые 2 и 3 сдвинуты к большим плотностям. Кривые 1–3 примерно сохраняют симметрию относительно плотности 0.5.

Рисунок 3 демонстрирует, что учет ФМ + КХП понижает критическую температуру, и это понижение тем больше и ближе к точному значению, чем больше размер фрагмента (кривые 1 на рис. 3а–в). С другой стороны, для любого размера фрагмента калибровка позволяет получить точное значение (βε)_c = 1.764, и для него можно рассматривать квазиточные значения термодинамических функций.

Наличие двух разных плотностей в системе (узлы фрагмента и узел КХП) вносит искажение в симметрию задачи – кривая расслаивания становится асимметричной и имеет критическую плотность отличную от значения θ_c = 0.5. При сохранении величин молекулярных параметров системы приближение Т_с к точному значению сопровождается смещением от θ_c = 0.5. Вопрос о целесообразности строгого сохранения симметрии определяется целью использования калибровочных функций. С увеличением размера фрагмента доля вклада от узла КХП в термодинамические функции уменьшается. Процедура калибровки позволяет за счет варьирования двух параметров получить и точное значение (βε)_c и θ_c = 0.5. Это (частично или практически полностью) нивелирует искажение симметрии купола фазовой диаграммы.

Корреляторы ФМ + КХП в надкритической области

На рис. 4 представлены изотермы корреляторов парного (рис. 4а) и пятого порядка (рис. 4б) в ФМ + КХП без калибровки при K = 4 и Q_g = 0, нормированных на аналогичные корреляторы в ФМ, ζ_n = θ_(n)(ФМ + КХП)/θ_(n)(ФМ), при τ = 2.5 (1, 5), 2.1 (2, 6), 1.7 (3, 7) и 1.3 (4, 8) для фрагментов 3 × 3 (1–4) и 4 × 4 (5–8).

Рис. 4.

Изотермы корреляторов парного (а) и 5-го (б) порядка для ФМ + КХП без калибровки.

Согласно рис. 4 корреляторы в ФМ + КХП выше ФМ и с ростом фрагмента они убывают, но это связано исключительно с завышенным координационным числом (к.ч.) для узлов фрагмента по сравнению с z на число K/N_fr, которое убывает с ростом размера фрагмента. Значение отношения корреляторов с ростом плотности убывает, достигая значения единицы, т.е. при полном заполнении системы мы получаем общий результат в ФМ + КХП, ФМ и КХП (отношение ФМ и КХП на рис. 1 также достигало значения единицы при высоких плотностях).

Чем выше порядок корреляторов, тем больше разница между ФМ + КХП и ФМ: на рис. 4а для ζ₂ максимальные значения составляют 1.17, а для ζ₅ на рис. 4б – 1.5. В пределе разреженного газа все корреляции отсутствуют и все методы дают единицу в точке θ = 0. Рисунок 4 показывает, как с ростом порядка корреляторов растет расхождение между ФМ и ФМ + КХП.

Аналогичные кривые с другой нормировкой были представлены на рис. 1, чтобы подчеркнуть отличия точных расчетов от кривых в КХП, а здесь демонстрируются отличия квазиточных величин от точных. Расчеты показывают, что в методе ФМ + КХП концентрационные зависимости всех корреляторов отражаются достаточно точно – отношения их величин находятся вблизи единицы по сравнению с точными значениями.

Квазиточные кривые расслаивания для кубической решетки

Предложенный подход к комбинированному использованию ФМ + КХП может быть применен к другим двух- и трехмерным решеткам. На рис. 5а показаны фазовые диаграммы для кубической решетки: точечной линией для КХП (1), пунктирными линиями для ФМ + КХП при K = 6, Q_g = 0 (без калибровки (2–4)) и сплошными линиями для ФМ + КХП с калибровкой (5–7) для фрагментов 2 × 2 × 2 (2, 5), 2 × 2 × 3 (3, 6) и 2 × 3 × × 3 (4, 7). В КХП (1) в критической точке имеем: (βε)_c = 0.812, θ_c = 0.5.

Рис. 5.

Фазовые диаграммы (а) и изотермы (б) трехмерной системы.

Так же, как и в двухмерной системе (кривая 3 на поле 2б), в трехмерной системе с увеличением размера фрагмента величина Т_c уменьшается, и все они находятся ниже точного значения ((βε)_c = = 0.888): фрагмент 2 × 2 × 2 (βε)_c = 1.019, θ_c = = 0.441, фрагмент 2 × 2 × 3 (βε)_c = 1.054, θ_c = 0.433, фрагмент 2 × 3 × 3 (βε)_c =1.060, θ_c = 0.428. Поэтому необходимо использовнаие калибровочных функций.

Кривые 5–7 отвечают калиброванным значениям K = 9.12 и Q_g = –0.74 для 2 × 2 × 2, K = 10.60 и Q_g = –0.74 для 2 × 2 × 3 и K = 11.74 и Q_g = –0.76 для 2 × 3 × 3 под точное значение (βε)_c = 1.764 и плотность θ_c = 0.5. Здесь с ростом размера фрагмента величина параметра K увеличивается, отдаляясь от K = z = 6. Величина Q_g  уменьшается с ростом размера фрагмента.

На рис. 5б показаны изотермы точечными линиями для КХП (1, 8), пунктирными линиями для ФМ (2, 5, 9, 12), тонкими сплошными линиями для ФМ + КХП без калибровки (3, 6, 10, 13) и толстыми сплошными линиями для ФМ + КХП с калибровкой (4, 7, 11, 14) для фрагментов 2 × 2 × 2 (2–4, 9–11) и 2 × 3 × 3 (5–7, 12–14) при температуре ниже критики τ = 0.7 (1–7) и выше 1.3 (8–14).

Рис. 5б демонстрирует, что изотермы КХП (1, 8) и ФМ (2, 5, 9, 12) при заданной температуре имеют общий центр в точке плотности θ = 0.5 и давления, отвечающего равновесному в КХП при расслаивании. Учет ФМ + КХП (3, 6, 10, 13) уменьшает давление, а учет корреляций (4, 7, 11, 14) еще больше уменьшает давление и уменьшает наклон изотерм в центре к оси плотности. Также давления уменьшаются с ростом фрагмента от 2 × × 2 × 2 (2–4, 9–11) к 2 × 3 × 3 (5–7, 12–14).

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

Точные значения термодинамических функций играют важную роль для тестирования методов моделирования. Они позволяют определить корректность и точность используемых приближенных подходов, что особенно важно для неидеальных систем. ФМ представляет собой прямой расчет статсуммы малой системы [12–14, 17]. Рассматриваемый фрагмент представляет собой периодический участок макроскопической поверхности. В нем отсутствуют какие-либо приближения, но в силу малости числа узлов системы он не в состоянии описывать фазовые переходы, которые являются следствием макроскопического числа частиц. Его внутренняя топография для неоднородных систем может быть произвольной.

Комбинированное использование ФМ и КХП в рамках рассмотренного варианта с формальным использованием узла КХП в качестве калибровочной функции существенно расширяет область точного описания физико-химических свойств в широких диапазонах температур и плотностей. ФМ + КХП позволяет использовать известные точные значения для Т_cr [5–9], для нахождения квазиточных функций распределений (унарных, парных и корреляторов более высоко порядка) в нужных диапазонах температур и плотностей, и получения квазиточных значений термодинамических функций. В частности, подход ФМ + КХП расширяет возможности калибровочных функций [15–17], которые ранее строились только на основе информации о величинах Т_c, на совместное использование величины Т_c и точных расчетов в надкритической области.

Возможность калибровки в ФМ + КХП фрагментов практически любого размера на точное значение Т_c  дает преимущества такого подхода по сравнению вариационным кластерным методом [20, 21] за счет ускорения расчетов (не менее чем на три порядка) и возможностью рассматривать большие по размеру фрагменты, что особенно важно для неоднородных систем.

Предложенный подход расширяет круг систем, для которых повышается точность расчета: это получение точных значений для фрагментов N_fr ≤ 30, возможность описания фазовых переходов до N_fr ≤ 24–25, а также возможность замены КХП в прикладных задачах с включением областей фрагментов порядка N_fr ≤ 10–16. Данные размеры определяют также размеры разных подсистем более общих задач, в которых возможно применение нового подхода ФМ + КХП для расчета или тестирования других численных подходов, что позволяет описывать физико-химические свойства многих реальных задач.

Работа выполнена при поддержке Российского фонда фундаментальных исследований (код проекта № 18-03-00030а).