Non-covalent interaction of carbon, silicon and germanium atoms

Full Text

1. ВВЕДЕНИЕ

Нековалентные взаимодействия проявляются в газах и жидкостях, а также при их взаимодействиях с поверхностями твердых тел [1]. Кроме того, нековалентные взаимодействия имеют решающее значение в образовании супрамолекулярных структур, которые встречаются в химии, биологии и нанотехнологиях [2]. Эти взаимодействия характеризуются малыми значениями энергии и времени образования ассоциатов, что обеспечивает возможность эффективного управления их структурой и физическими свойствами.

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

К настоящему времени накоплен достаточно обширный материал по эмпирическим парным потенциалам [11]. Параметры этих потенциалов оцениваются по вириальным коэффициентам, температурам и теплотам фазовых переходов, модулям упругости, коэффициентам диффузии, вязкости, теплопроводности и теплового расширения, а также по данным измерений такими методами, как атомно-силовая спектроскопия, аппарат поверхностных сил и спектроскопия ядерного магнитного резонанса. Выбор того или иного способа оценки параметров диктуется спецификой решаемой задачи. При этом численные значения параметров парных потенциалов для одних и тех же систем могут существенно различаться.

Среди теоретических ab initio методов расчета потенциалов взаимодействия на сегодняшний день центральную позицию занимают методы, основанные на теории функционала электронной плотности [12]. К ним относятся такие приближения, как Томаса–Ферми (Thomas–Fermi( TF)), Томаса–Ферми–Дирака (Thomas–Fermi–Dirac (TFD)), Гордона–Кима (Gordon-Kim (GK)), Кона–Шэма (Kohn–Sham (KS)), Пердью–Бурка–Эрнзерхофа (Perdew–Burke–Ernzerhof (PBE)). Широко также применяется модифицированный метод погруженного атома (modified embedded-atom method (MEAM)). Основная задача всех этих методов заключается в получении надежно охарактеризованных межатомных потенциалов.

Расчетам потенциалов ковалентного взаимодействия для элементов подгруппы углерода (C, Si, Ge) посвящено большое количество работ, среди которых можно отметить работы [13–15]. Нековалентное взаимодействие, представляющее собой дисперсионное взаимодействие Ван-дер-Ваальса, для данной подгруппы в рамках метода функционала электронной плотности изучалось в работах [16, 17]. Были определены константы взаимодействия C₆ только для взаимодействия одинаковых атомов. В работе [17] расчеты проведены как для свободных атомов, так и для атомов в составе твердого вещества. В отличие от [16, 17], целями настоящей статьи были: 1) расчет из первых принципов потенциалов нековалентного взаимодействия свободных атомов углерода, кремния и германия во всевозможных сочетаниях; 2) вычисление параметров потенциалов Леннарда-Джонса и Морзе; 3) проверка эмпирических правил комбинирования параметров потенциалов Лоренца–Бертло; 4) расчет констант дисперсионного взаимодействия и вириальных коэффициентов. Хорошо известно, что действие дисперсионных сил притяжения предшествует ковалентному связыванию атомов. По этой причине в реальных системах на потенциальных кривых могут наблюдаться несколько минимумов, которые отвечают нековалентному и ковалентному типам взаимодействий. Основной результат статьи заключается в расчете именно дисперсионного взаимодействия, которое имеет место при больших расстояниях между атомами, т.е. когда это расстояние, во всяком случае, превышает сумму ковалентных радиусов атомов. В расчетах используется метод функционала электронной плотности, а выражения для электронных плотностей берутся в приближении Хартри–Фока с учетом оболочечной структуры взаимодействующих атомов.

2. ТЕОРЕТИЧЕСКАЯ МОДЕЛЬ

В данном разделе во всех формулах используется атомная система единиц ($\hslash =e=m_e=1$). В качестве теоретической модели для нековалентного взаимодействия двух атомов примем приближение Гордона–Кима [18, 19]. В рамках приближения предполагается, что результирующая электронная плотность равна алгебраической сумме плотностей отдельных атомов. Кулоновская составляющая потенциала взаимодействия включает в себя отталкивание между ядрами, отталкивание между электронными облаками и притяжение между ядрами и электронными облаками атомов. Все эти части группируются в одно выражение:

$U_{\text{I}}\left(r\right)=16{\pi }^2\underset{0}{\overset{\infty }{}}\underset{0}{\overset{\infty }{}}{\rho }_a\left(x\right){\rho }_b\left(y\right)I\left(r,x,y\right)x^2{y}^{2}dxdy,$ (1)

$\begin{array}{c}I\left(r,x,y\right)=\frac{1}{r}+F\left(r,x,y\right)-\\ -\frac{2}{r+x+\left|r-x\right|}-\frac{2}{r+y+\left|r-y\right|},\end{array}$ (2)

$F\left(r,x,y\right)=\left\{\begin{array}{l}\frac{2}{r+x+\left|r-x\right|}, y<\left|r-x\right|,\\ \frac{1}{2}\left(\frac{1}{x}+\frac{1}{y}\right)-\frac{r}{4xy}-\frac{{\left(x-y\right)}^2}{4rxy}, \\ \left|r-x\right|<y<\left|r+x\right|,\\ \frac{1}{y}, y>r+x,\end{array}\right.$ (3)

где r – расстояние между атомами, $\rho$_a,b – электронные плотности атомов, индексы “a” и “b” обозначают взаимодействующие атомы. Квантовостатистическая составляющая, относящаяся к электронному газу, вычисляется с помощью интеграла

$\begin{array}{c}{U}_{II}\left(r\right)=\frac{\pi r^3}{4}\underset{1}{\overset{\infty }{}}\underset{-1}{\overset{1}{}}\left[\left({\rho }_a+{\rho }_b\right)E\left({\rho }_a+{\rho }_b\right)-\right.\\ -\left.{\rho }_{a}E\left({\rho }_a\right)-{\rho }_{b}E\left({\rho }_b\right)\right]d\mu d\lambda ,\end{array}$ (4)

$E\left(\rho \right)=\frac{3}{10}{\left(3{\pi }^2\rho \right)}^{2/3}-\frac{3}{4}{\left(\frac{3}{\pi }\rho \right)}^{1/3}+E_c\left(\rho \right) ,$ (5)

$E_c\left(\rho \right)=\left\{\begin{array}{l}2.711816{\rho }^{1/2}-0.706052{\rho }^{1/3}-\\ -3.819822{\rho }^{2/3}-1.319676{\rho }^{5/6},\rho <{\eta }_1,\\ -0.7062-0.00633\ln \left(\rho \right), {\eta }_1\le \rho \le {\eta }_2,\\ -0.010367\ln \left(\rho \right)-\left(0.001861\ln \left(\rho \right)+\right.\\ +\left.0.008869\right){\rho }^{-1/3}-\mathrm{0.062849,}\rho >{\eta }_2,\end{array}\right.$

$\begin{array}{l}{\rho }_a\left(r;\lambda ,\mu \right)={\rho }_a\left[\frac{r\left(\lambda +\mu \right)}{2}\right], \\ {\rho }_b\left(r;\lambda ,\mu \right)={\rho }_b\left[\frac{r\left(\lambda -\mu \right)}{2}\right],\end{array}$

где отдельные слагаемые в (5) отвечают кинетическому, обменному и корреляционному вкладам, соответственно; $\eta$₁ = 2.387325 · 10^-4 и $\eta$₂ = 0.696013 – параметры сшивки в корреляционном вкладе. Двойной интеграл в (4) выражен через эллипсоидальные координаты ($\lambda , \mu$). При его вычислении с учетом (5) под p понимаются электронные плотности взаимодействующих атомов p_a,b и их сумма p_a + p_b. В настоящей работе было найдено, что для корреляционной энергии вполне применима полуэмпирическая (второй сомножитель воспроизводит правильную асимптотику по p^-1/3 [18, 19]) функция вида

$\begin{array}{c}{E}_c\left(r_s\right)=\left[1-\text{exp}\left(-0.0662r_s^{2.123}\right)\right]\times \\ \times \left(-\frac{0.438}{r_s}+\frac{1.325}{r_s^{1.5}}-\frac{1.47}{r_s^2}-\frac{0.4}{r_s^{2.5}}\right), \end{array}$

$r_s={\left(\frac{3}{4\pi \rho }\right)}^{1/3},$

где r_s – радиус Вигнера–Зейтца. По отношению к выражению (20) из работы [18] ошибка данной аппроксимации при 0.01 ≤ r_s ≤ 10 не превышает 5%, а при r_s > 10 она фактически равна нулю. Результирующая потенциальная энергия взаимодействия двух атомов есть

$U\left(r\right)=U_{\text{I}}\left(r\right)+U_{\text{II}}\left(r\right).$ (6)

В модели (1)–(6) требуется задать p_a,b. Наиболее корректно с физической точки зрения это можно сделать [20], если использовать в расчетах квантовомеханические электронные плотности атомов, вычисленные методом Хартри–Фока. В работе [21] получена формула

$\begin{array}{c}\rho \left(\xi \right)=\frac{Z}{4\pi }\times \\ \times \left[\underset{i}{}{A}_i\exp \left(-{\alpha }_i\xi \right)+\underset{j}{}{B}_j\left({\beta }_j\xi -2\right)\exp \left(-{\beta }_j\xi \right)\right],\end{array}$ (7)

где Z – заряд ядра атома; A_i, a_i, B_j и b_j – числовые коэффициенты, которые определяются по данным из [21]; x – координата, отсчитываемая от ядра атома. Отметим, что формула (7) успешно применялась для расчета диамагнитной восприимчивости и атомных факторов рассеяния рентгеновского излучения.

Численные расчеты по формулам (1)–(7) проводились следующим образом. Обычно для приближенного вычисления несобственных интегралов с одним бесконечным пределом применяется метод Гаусса–Лагерра [22, стр. 890]. В нашем случае подынтегральные функции в формулах (1) и (4) имеют сложный осциллирующий характер поведения на коротких расстояниях, вследствие чего достижение приемлемой точности требует вычисления порядка ста и более корней полиномов Лагерра. Вдобавок большая часть из них приходится на участки, где значения функций ничтожно малы и поэтому не привносят ощутимого вклада в результирующую сумму. В связи с этим метод Гаусса–Лагерра оказывается не столь результативным, несмотря на рекомендацию из работы [18]. Нами интегрирование проводилось в ограниченных областях, но в пределах которых располагаются “эффективные” части функций. Область интегрирования для (1) представляла собой квадрат размером 60×60, а для (4) – прямоугольник размером 2×60. При этом расстояния в 40–60 а.е. уже достаточно велики, чтобы служить хорошим приближением к бесконечности. Расстояние между атомами во всех расчетах изменялось в широком интервале: от 0.5 до 20 а.е. с шагом 0.05 а.е.

Численное интегрирование проводилось методом трапеций [22, стр. 885], который наиболее устойчив к накоплению локальных ошибок. Расчетная сетка состояла из ≈5 · 10⁸ точек, что в итоге гарантировало 4–6 верных знаков после запятой. Сам алгоритм был распараллелен для проведения вычислений на графическом процессоре и реализован в программной среде CUDA C++. Это позволило уменьшить время работы программы на порядки. Разработанные численный алгоритм и программа тестировались путем сравнения расчетов потенциала межатомного взаимодействия с данными из работ [23, 24] для системы Ar–Ar (рис. 1).

 

Рис. 1. Потенциал взаимодействия для системы Ar–Ar: сплошная кривая – настоящая работа, штриховая – работа [23], точечная – работа [24].

 

Параметры модельных парных потенциалов определялись по табулированным кривым зависимости U(r). Минимумы потенциальных функций находились напрямую из табличных данных посредством их перебора. При этом погрешность по расстоянию соответствует шагу разбиения и в нашем случае не превышает 0.05 а.е. Показатели экспонент также оценивались численно по графикам функций U(r) в точках минимума с приближенным согласованием по жесткости межатомной связи k_a, т.е. через вторые производные функций в точках минимума. Данный способ оказался более эффективным по сравнению с такими опробованными нами оптимизационными методами, как методы Ньютона, Левенберга–Марквардта и сопряженных градиентов, реализация которых критически зависит от начальных условий и не всегда обеспечивает единственное решение.

3. АНАЛИЗ РЕЗУЛЬТАТОВ

Результаты проведенных расчетов в виде графиков представлены на рис. 2–4. Расчеты позволяют идентифицировать параметры модельных парных потенциалов, а также проверить правила комбинирования Лоренца–Бертло, которые широко используются на практике [5]. Параметры этих потенциалов всегда определяются ad hoc и зависят от особенностей той или иной задачи (см., например, [25–28]). В простейшем случае для парного межатомного потенциала взаимодействия имеем аппроксимацию [10]

$U\left(r\right)=D\left(X^2-2X\right),$ (8)

где D – глубина потенциальной ямы, X = X(r) – функция, которая задает зависимость от расстояния между атомами. Если принять X = (r₀/r)⁶ и X = exp{a(r₀ – r)}, то из (8) получим соответственно потенциалы Леннарда–Джонса и Морзе:

$U\left(r\right)=D\left[{\left(\frac{r_0}{r}\right)}^{12}-2{\left(\frac{r_0}{r}\right)}^6\right],$ (9)

$k_a=U^{\text{'}\text{'}}\left(r_0\right)=\frac{72D}{r_0^2},$

$U\left(r\right)=D\left[\exp \left\{2\alpha \left(r_0-r\right)\right\}-2\exp \left\{\alpha \left(r_0-r\right)\right\}\right],$ (10)

$k_a=U^{\text{'}\text{'}}\left(r_0\right)=2D{\alpha }^2,$

где r₀ – равновесное расстояние между атомами, a определяет кривизну потенциальной функции в области минимума. Потенциал (9) при r → ∞ имеет физически правильную степенную асимптотику, свойственную силам Ван-дер-Ваальса. Однако наши расчеты показывают, что расчет по формуле (9) дает слишком крутой рост энергии отталкивания при уменьшении межатомного расстояния. В этом отношении более оправданным оказывается потенциал (10), в котором учитывается плавность силы отталкивания за счет параметра a. Но недостатком формулы (10) является отсутствие особенности при r → 0. В определенной мере положительные свойства потенциалов (9) и (10) можно объединить, если в (8) принять функцию вида X = (r₀/r)^mexp{a(r₀ – r)}, где m ≥ 0. В итоге находим

$\begin{array}{c}U\left(r\right)=D\left[{\left(\frac{r_0}{r}\right)}^{2m}\exp \left\{2\alpha \left(r_0-r\right)\right\}-\right.\\ -\left.2{\left(\frac{r_0}{r}\right)}^m\exp \left\{\alpha \left(r_0-r\right)\right\}\right],\end{array}$ (11)

$k_a=U^{\text{'}\text{'}}\left(r_0\right)=\frac{2D}{r_0^2}{\left(\alpha r_0+m\right)}^2.$

При m = 0 из (11) следует формула (10) для потенциала Морзе. Экранированный кулоновский потенциал для энергии отталкивания отвечает значению m = 1/2. При m = 1 из (11) получается потенциал, предложенный ранее одним из авторов для моделирования силовых сигналов в атомно-силовом микроскопе в работе [28]; в области межатомного притяжения здесь возникает потенциал типа потенциала Юкавы. Таким образом, формула (11) более гибко описывает расчеты в приближении электронного газа во всем диапазоне расстояний (при r ∈ (0,∞)), чем формулы (9) и (10) и может быть полезной при исследовании свойств реальных газов и жидкостей методами молекулярной динамики и Монте-Карло. Кроме того, потенциал (11) представляется перспективным для создания упрощенных численных моделей адсорбционных явлений в наноразмерных системах, которые эффективно изучаются с применением теории функционала электронной плотности [29–33].

В табл. 1 приведены параметры потенциалов (9)–(11), вычисленные по (1)–(7) в заданном интервале межатомных расстояний. Обращает на себя внимание то, что параметр a в выражениях (10), (11) слабо изменяется и имеет приблизительно одно и то же значение для различных атомов – соответственно 1.4 и 1.2 Å^-1. Для сравнения приведены также суммы ковалентных радиусов r_c для взаимодействующих атомов [34]. Их значения задают расстояния, на которых образуются ковалентные неполярные связи. При таких расстояниях (~0.2 нм) для атомов подгруппы углерода использованное приближение уже не работает.

 

Таблица 1. Параметры модельных потенциалов

Атомная пара	D, эВ	r0, Å	rc, Å	a, Å-1		m
				(10)	(11)
С–С	0.021	3.286	1.50	1.473	1.219	1/2
Si–Si	0.088	3.250	2.32	1.458	1.204	1/2
Ge–Ge	0.071	3.685	2.42	1.402	1.165	1/2
C–Si	0.035	3.467	1.91	1.394	1.148	1/2
C–Ge	0.035	3.540	1.96	1.45	1.206	1/2
Si–Ge	0.078	3.491	2.37	1.413	1.166	1/2

 

Наиболее часто параметры взаимодействия разнородных атомов оцениваются по параметрам взаимодействия одинаковых атомов согласно правилам комбинирования Лоренца–Бертло [5]:

$D_{ab}={\left(D_a{D}_b\right)}^{1/2}, r_{ab}=\frac{r_a+r_b}{2}.$ (12)

Из теории функционала электронной плотности следует, что правила комбинирования (12) не могут считаться универсальными. Так, в нашем случае для системы C–Si эти правила не выполняются даже качественно. Расчеты из первых принципов дают r_ab > r_a,r_b, что не согласуется с (12). Это объясняется электронным строением взаимодействующих атомов. Вдали от ядер имеется область пространства, в которой электронная плотность атома углерода превышает электронную плотность атома кремния (см. рис.1), что приводит к изменению энергии отталкивания.

Расчеты потенциальных функций по формулам (1)–(7) в области действия сил притяжения позволяют провести оценки дисперсионных констант C₆ с помощью метода наименьших квадратов и функции U(r) = –C₆/r ⁶. С другой стороны, имеет место формула Слэтера–Кирквуда, которая получена вариационным методом [1, стр.36]:

$C_6=\frac{3}{2}\frac{{\alpha }_a{\alpha }_b}{{({\alpha }_a/N_a)}^{1/2}+{({\alpha }_b/N_b)}^{1/2}},$ (13)

где a_a,b и N_a,b – статические поляризуемости и валентности атомов. В нашем случае были приняты следующие значения: a_C = 12 а.е., a_Si = 34.6 а.е., a_Ge = 41.7 а.е. [16] и N_a = N_b = 4. В табл. 2 приведены результаты наших расчетов в сравнении с данными из работ [16, 17]. Можно видеть, что все значения согласуются друг с другом по порядку величины.

 

Таблица 2. Константы дисперсионного взаимодействия

Атомная пара	С6, а.е.
	(1)–(7)	(13)	[16]	[16, 17]
C–C	42.3	62.4	42.7	46.6
Si–Si	340	305	317	305
Ge–Ge	452	404	359	354
C–Si	128	133	–	–
C–Ge	147	151	–	–
Si–Ge	394	351	–	–

 

Полученными результатами можно воспользоваться для расчета второго вириального коэффициента [35, 36]:

$B\left(T\right)=-2\pi {\int }_0^{\infty }\left[\text{exp}\left(-\frac{U\left(r\right)}{k_{B}T}\right)-1\right]r^{2}dr,$ (14)

где T – температура, k_B – постоянная Больцмана, потенциал U(r) определяется по (1)–(7) или (9)–(11). C физико-химической точки зрения данный расчет не столь информативен по отношению к газообразным фазам углерода, кремния и германия в широком температурном диапазоне, но позволяет осуществить выбор модельного потенциала применительно к вириальному уравнению состояния, т.е. к моделированию термодинамических свойств газов Ван-дер-Ваальса.

 

Рис. 2. Распределение радиальных электронных плотностей атомов.

 

Рис. 3. Потенциалы взаимодействия гомоатомных пар.

 

Рис. 4. Потенциалы взаимодействия гетероатомных пар.

 

Реперный расчет функции B(T) для системы Ar–Ar демонстрирует удовлетворительное совпадение с литературными данными. Для глубины потенциальной ямы и равновесного расстояния в нашем случае имеем D/k_B = 129.51 K и r₀ = 3.649 Å. В [18] были найдены значения D/k_B = 126.752 и r₀ = 3.63 Å. Некоторое несовпадение между нашими расчетами и расчетами из [18] по системе Ar–Ar обусловлено выбором функции электронной плотности. Оптимизация согласно экспериментальным данным и формуле (8) из [36], которая получается путем интегрирования (14) с учетом потенциала Леннарда-Джонса (9), дает значения того же порядка: D/k_B = 119 К и r₀ = 3.859 Å. Разница между нашими расчетами функции B(T) и расчетами по формуле (8) из [36] в небольшой окрестности точки Бойля и ниже этой точки не превышает 10%. Расчеты с потенциалом из [24] с экспериментальными данными для аргона согласуются хуже (см. [37]).

На рис. 5 представлены результаты расчета функции B(T) для газа Ван-дер-Ваальса атомов углерода. Во всех случаях интеграл (14) рассчитывался в интервале, в котором табулировалась функция U(r) в приближении электронного газа (1)–(7). Численное вычисление интеграла (14) проводилось методом трапеций. Как и ожидалось, наиболее близкими к модели электронного газа оказываются потенциалы (10) и (11). Если же варьировать параметры модельных потенциалов (9)–(11), то во всех случаях можно добиться лучшего совпадения с расчетами из первых принципов или экспериментальными данными по второму вириальному коэффициенту, как это сделано, например, в работе [36]. Но тогда сами потенциальные функции будут заметно отличаться друг от друга и от фундаментальной зависимости, которая строится по формулам (1)–(7).

 

Рис. 5. Температурная зависимость второго вириального коэффициента для углерода.

 

4. ЗАКЛЮЧЕНИЕ

В настоящей работе построены зависимости потенциальной энергии нековалентного взаимодействия от расстояния для атомов углерода, кремния и германия во всех сочетаниях. В известных нам монографиях и статьях, посвященных этим и другим близким вопросам, подобные зависимости отсутствуют. Использовалась теория электронного газа [18–20], в которой электронная плотность атомов учитывалась в приближении Хартри–Фока [21]. Определены параметры модельных межатомных потенциалов, которые могут использоваться в различных статистических и термодинамических теориях. Показано, что широко используемые правила комбинирования Лоренца–Бертло для параметров взаимодействия не всегда работают. Правила, в частности, дают ошибочные результаты для системы C–Si. В качестве примера проведены модельные расчеты второго вириального коэффициента для пара атомов углерода. Эти расчеты, разумеется, нельзя считать физически приемлемыми в широком температурном диапазоне, особенно в области низких температур, соответствующих конденсации пара, но они формально подтверждают применимость с почти одинаковым успехом всех проанализированных модельных потенциалов для построения уравнений состояния.

Авторы искренне признательны Рецензенту за высказанные замечания, устранение которых позволило улучшить статью.

About the authors

A. A. Sokurov

Institute of Applied Mathematics and Automation – the filial branch of Federal State Budgetary Scientific Establishment “Federal Scientific Center “Kabardin-Balkar Scientific Center of Russian Academy of Sciences”

Email: rsergo@mail.ru
Russian Federation, Nalchik

S. S. Rekhviashvili

Institute of Applied Mathematics and Automation – the filial branch of Federal State Budgetary Scientific Establishment “Federal Scientific Center “Kabardin-Balkar Scientific Center of Russian Academy of Sciences”

Author for correspondence.
Email: rsergo@mail.ru
Russian Federation, Nalchik

References

Supplementary files