ЖЭТФ, 2022, том 161, вып. 1, стр. 86-103

МЕЖАТОМНОЕ ВЗАИМОДЕЙСТВИЕ НА ГРАНИЦЕ

АЛЮМИНИЙ-ФУЛЛЕРЕН C₆₀

В. В. Решетнякa*, О. Б. Решетняк^a, А. В. Аборкин^b, А. В. Филиппов^a,c

^a ГНЦ РФ Троицкий институт инновационных и термоядерных исследований

108840, Троицк, Москва, Россия

^b Владимирский государственный университет им. А. Г. и Н. Г. Столетовых

600000, Владимир, Россия

^c Объединенный институт высоких температур Российской академии наук

125412, Москва, Россия

Поступила в редакцию 28 июля 2021 г.,

после переработки 17 сентября 2021 г.

Принята к публикации 17 сентября 2021 г.

Предложена модель взаимодействия атомов алюминия и углерода на границе Al/C60. В рамках теории

функционала плотности рассчитаны энергия связи и положение фуллерена на подложке Al(111). Полу-

ченные результаты использованы для определения параметров потенциала Леннарда-Джонса, который

затем применялся в исследованиях методом молекулярной динамики. Проведенное теоретическое иссле-

дование десорбции фуллеренов с алюминиевой подложки показало хорошее согласие результатов с из-

вестными из литературы экспериментальными данными. Изучены капиллярные эффекты, возникающие

на границе между расплавом алюминия и погруженными в него фуллеренами. Положительное значение

удельной свободной энергии, приходящейся на единицу поверхности Al/C60, указывает на слабую сма-

чиваемость молекул расплавом. Рассчитанное значение времени диффузионной релаксации оказалось

примерно на два порядка величины меньше характерного времени коагуляции фуллеренов, что указывает

на наличие сил отталкивания между ними. Обсуждается активационный характер процесса коагуляции

и капиллярная природа взаимодействия между фуллеренами.

DOI: 10.31857/S0044451022010084

ве данных сканирующей туннельной микроскопии,

фотоэлектронной спектроскопии и абсорбционной

1. ВВЕДЕНИЕ

рентгеновской спектроскопии в этих работах были

предложены две различные модели реконструкции

Углеродные частицы различной структуры ши-

√

поверхности: 2

3×2

3 и 6 × 6. Установлено, что

роко используются при разработке композиционных

энергии связи различны для этих типов реконструк-

материалов с требуемыми свойствами [1-4]. В част-

√

ции. Структура слоя фуллеренов 2

3 × 2

3 ста-

ности, значительный интерес представляют алюмо-

бильна при более низких температурах, a при 490 K

матричные композиты, наполненные фуллеренами

покрытие приобретает структуру 6 × 6. Десорбция

C₆₀ [5-7].

оставшихся фуллеренов с поверхности происходит

Одним из ключевых вопросов при создании ком-

при 730 К.

позиционных материалов является обеспечение вза-

Согласно данным микроскопии [8] в реконструк-

имодействия матрицы и наполнителя. Взаимодей-

√

ции 2

3×2

3 существуют два различных способа

ствие фуллеренов C₆₀ с различными кристалличе-

расположения фуллеренов на поверхности Al(111).

скими поверхностями алюминия ранее изучалось

Для них характерны различные энергии связи и рас-

экспериментально [8-13] и теоретически [14]. Из ра-

стояния между молекулой и подложкой. Модель из

бот [8-11] известно, что при взаимодействии C₆₀

[14] объясняет этот эффект смещением молекулы

с поверхностью алюминия между атомами Al и C

в сторону подложки с вакансией. Последующие ис-

возникает ковалентная химическая связь. На осно-

следования, проведенные для поверхностей различ-

^* E-mail: viktor.reshetnyak84@gmail.com

ных металлов с плотной упаковкой атомов, показа-

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

ли, что фуллерены склонны располагаться на этих

рассчитана как аналитически, так и численно, с

поверхностях над вакансиями. При этом кристал-

применением метода молекулярной динамики (МД).

лическая структура металла определяет тип рекон-

Для ускорения процесса десорбции при этом бы-

струкции [15].

ли использованы заниженные значения энергии свя-

Взаимодействие между фуллеренами и поверх-

зи. Для лабораторных временных интервалов, со-

ностями Al(100) и Al(110) менее изучены. В статье

поставимых с характерными временами из экспе-

[9] показано, что при взаимодействии C₆₀ с Al(110)

риментов [13], применялась аналитическая модель

атомы углерода и алюминия образуют ковалент-

и значение энергии связи, полученное в настоя-

ные связи, а десорбция фуллеренов происходит при

щей работе в рамках теории функционала плот-

730 K. Авторы не наблюдали различных способов

ности (density functional theory, DFT). Сравнение

расположения молекул на этих поверхностях и со-

показывает, что предложенная аналитическая мо-

ответствующих пиков десорбции.

дель и потенциал межатомного взаимодействия вос-

В данной работе для описания взаимодействия

производят известные экспериментальные данные с

атомов Al и C для систем Al-C₆₀ использовался по-

высокой точностью. Результаты моделирования ме-

тенциал Леннарда-Джонса:

тодом МД достаточно хорошо согласуются с ана-

)₁₂

)₆]

литическими. Полученные данные обсуждаются в

[( σ

⁽σ

u(r) = 4ε

(1)

разд. 3.2.

В статьях [17-20] указывалось, что значения

энергии связи фуллеренов с подложкой, структу-

где ε — глубина потенциальной ямы, а значение па-

ры и температуры десорбции весьма близки для по-

раметра σ определяет положение нуля потенциаль-

верхностей некоторых различных металлов. В част-

ной кривой. Параметры потенциала вычислялись с

ности, полученные нами результаты для алюмини-

использованием аналитической модели взаимодей-

евой подложки очень близки к данным для подло-

ствия однородной сферы с подложкой и расчетов

жек из благородных металлов. Это позволяет наде-

ab initio. Подробное обсуждение предложенной мо-

яться, что обобщение результатов настоящей работы

дели приведено в разд. 3.1.

на случай систем Au-C₆₀ и Ag-C₆₀ не потребует су-

Для сравнения с экспериментом использовалась

щественной доработки. Однако возможность такого

аналитическая модель процесса температурно про-

обобщения требует дополнительного исследования,

граммируемой десорбции (ТПД), основанная на

которое выходит за рамки данной статьи.

приближениях теории переходного состояния. Неиз-

В настоящей работе также изучены капилляр-

вестные параметры модели определялись по предва-

ные явления на границах между расплавленным

рительно рассчитанным ab initio значениям рассто-

алюминием и погруженными в него фуллеренами.

яния и энергии связи молекулы C₆₀ с поверхностью

Выполнен расчет скорости коагуляции, результа-

Al(111).

ты которого указывают на наличие взаимодействия

В методе ТПД температура системы поднима-

между фуллеренами на больших расстояниях. Уста-

ется по определенному закону (обычно линейному)

новлено, что процессы коагуляции и диффузии фул-

в заданном диапазоне значений. Во время нагрева

леренов имеют разные характерные времена релак-

контролируется поток молекул с поверхности. Пре-

сации и могут быть разделены на «медленные» и

имуществом подхода является возможность деталь-

«быстрые». Давление Лапласа p и поверхностное на-

ного анализа десорбции, что позволяет с высокой

тяжение γ были рассчитаны по данным МД. Затем,

точностью определять энергии связи для различных

учитывая изменение объема фуллеренов под дей-

позиций молекул на подложке. Подробное описание

ствием капиллярных сил, был вычислен объемный

метода можно найти, например, в обзорной статье

модуль. Обсуждению результатов исследования ка-

[16]. В статье [13] метод ТПД был применен для

пиллярных явлений посвящен разд. 3.3.

исследования взаимодействия между фуллеренами

C₆₀ и кристаллическими алюминиевыми поверхно-

2. ОПИСАНИЕ МОДЕЛЕЙ И МЕТОДОВ

стями. Была определена зависимость плотности по-

ЧИСЛЕННОГО ИССЛЕДОВАНИЯ

крытия монослоя фуллеренов n от температуры T

при постоянной скорости нагрева 3 К/с.

2.1. Расчет ab initio взаимодействия между

В данной работе выполнено многомасштабное

фуллереном и поверхностью Al(111)

исследование кинетики десорбции в различных тем-

пературных диапазонах. Для временных интерва-

Расчеты в рамках DFT проводились с использо-

лов длительностью порядка 1 нс зависимость n(T )

ванием программы CP2K [21,22]. Расчет энергии и

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

Рассматривались два варианта расположения

фуллерена на алюминиевой подложке. В первом

случае молекула C₆₀ размещалась над одним из ато-

мов алюминия, согласно модели hollow-6 из статьи

[14]. Если на подложке нет вакансии, то это положе-

ние соответствует наиболее стабильному состоянию

системы. Во втором случае фуллерен располагался

над вакансией на поверхности алюминия. Относи-

тельное положение молекулы C₆₀ и атомов алюми-

ния схематично показано на рис. 1.

2.2. Атомистическое моделирование с

использованием эмпирических силовых

полей

Все расчеты с эмпирическими межатомными по-

тенциалами были выполнены с использованием про-

граммы LAMMPS [33, 34]. Для учета взаимодей-

Рис. 1. Схема расположения атомов при взаимодействии

ствия между атомами углерода применялся потен-

фуллерена C60 с поверхностью Al(111). а — конфигурация

hollow-6, подложка без вакансии; б — фуллерен располо-

циал Терсоффа [35,36]. Для атомов алюминия была

жен над вакансией на подложке

использована модель погруженного атома (EAM) с

параметризацией, предложенной в [37]. Потенциал

Леннарда-Джонса вида (1) использовался для опи-

сания взаимодействия между атомами углерода и

сил выполнялся с учетом всех электронов в рамках

алюминия. Расчет параметров потенциала обсужда-

метода гауссовых и присоединенных плоских волн

ется в разд. 3.1. Интегрирование уравнений движе-

[23] с использованием базисного набора pob-TZVP

ния осуществлялось численно с шагом dt = 0.2 фс.

[24]. Обмен и корреляция электронов учитывались

При расчете десорбции фуллеренов с поверх-

в модели Пердью - Бурке - Эрнцерхофа [25].

ности алюминия были рассмотрены кристалличес-

Позиции атомов и параметры гранецентриро-

кая Al(111) и расплавленная алюминиевые пленки.

ванной кубической (ГЦК) кристаллической ячей-

Пленка толщиной в 10 атомных слоев была помеще-

ки алюминия были предварительно оптимизирова-

на параллельно плоскости xy. Размеры пленки опре-

ны в рамках DFT. Расчетное значение параметра

делялись суперячейкой 36×36 для кристаллической

решетки алюминия a = 4.0775Å совпало с экспе-

подложки Al(111). На границах вычислительной об-

риментальным 4.0469Å [26] в пределах 1 %. Полу-

ласти вдоль осей x и y были использованы периоди-

ченные данные далее использовались для построе-

ческие граничные условия (ГУ). На верхнем крае

ния пленки Al(111) с орторомбической суперячей-

z = z_max, расположенном на расстоянии 10 нм над

кой размером 5 × 6 и толщиной в четыре атомных

верхним атомным слоем, задавалось ГУ типа «от-

слоя. Пленка располагалась в плоскости xy, рассто-

ражающая стенка». Взаимодействие атомов с ниж-

яния между поверхностными атомами и граница-

ней границей определялось потенциалом Леннар-

ми расчетной ячейки в направлении z были зада-

да-Джонса 9-3 вида

ны равными 10Å. Интегралы по зоне Бриллюэна

[

]

)₉

)₃

вычислялись с использованием сетки Монхорста -

2 (σ

⁽σ

πn₁εσ³

(2)

Пака 3 × 3 × 1 [27]. Для предварительного тести-

рования модели были рассчитаны значения энер-

где n₁ — средняя плотность атомов алюминия в под-

гии поверхности E_s и энергии образования вакан-

ложке. В качестве параметров для алюминия бы-

сии E_v. Значение E_s = 1.06 эВ хорошо согласуется

ли взяты следующие значения: ε = 0.174 эВ, σ =

с результатами расчетов ab initio [28-30], где бы-

= 2.925Å [38]. Радиус обрезки потенциала выбирал-

ли получены значения в диапазоне от 0.67 эВ до

ся равным r_c = 2.5σ ≈ 7.3Å. Поскольку толщина

1.46 эВ, и экспериментальными данными из статьи

пленки в расчетах составляла около 25Å, действие

[31], E_s = 1.14-1.18 эВ. Рассчитанная энергия обра-

потенциала ограничивалось несколькими ближай-

зования вакансии E_v = 0.65 эВ согласуется с данны-

шими слоями атомов алюминия, а атомы углерода

ми [32], где было получено значение 0.61 эВ.

не взаимодействовали с нижней границей ячейки.

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

Можно получить потенциал в виде (2), рассчитав

энергию взаимодействия между атомом и полубес-

конечной однородной подложкой, предполагая при

этом, что потенциал межатомного взаимодействия

имеет вид (1) [39]. Для этого необходимо вычислить

следующий интеграл по объему подложки:

∞

∫

w = 2πn₁

u(L)L(L - H) dL.

(3)

Здесь L — расстояние от атома до точек элемента

объема dV (L), полученного разбиением подложки

сферическими поверхностями радиусами L и L+dL,

как показано на рис. 2a.

Моделирование ТПД проводилось с поверхно-

стями алюминия в кристаллическом и жидком со-

стояниях. Пленка алюминия предварительно стаби-

лизировалась при заданной начальной температуре

400 K в первом случае и 1000 K — во втором. Затем

на поверхность помещались 108 фуллеренов в соот-

√

ветствии с моделью реконструкции 2

3×2

3, по-

сле чего система стабилизировалась в течение 20 пс.

После стабилизации система нагревалась до целе-

вой температуры 800 K для кристалла и 2200 K для

жидкости. Продолжительность нагрева составляла

1 нс. Контроль температуры осуществлялся с ис-

пользованием термостата Ланжевена с параметром

температурной релаксации 20 фс. Для обеспечения

десорбции в заданных временных масштабах и тем-

пературном диапазоне значение параметра связи ε

потенциала (1) для атомов углерода и алюминия бы-

ло занижено, моделирование десорбции с поверхно-

сти кристаллической подложки было выполнено для

значения ε₁, уменьшенного в четыре раза относи-

тельно исходного ε, а для десорбции с поверхности

расплава было задано значение ε₂ = 0.7ε. Всего бы-

ло выполнено пять статистически независимых рас-

четов, по результатам которых вычислялась усред-

ненная зависимость плотности покрытия от темпе-

ратуры n(T ).

Для расчета спектральной плотности колебаний

центров масс молекул была рассмотрена система из

174960 атомов алюминия и углерода. Система вклю-

чала пленку Al(111) из 108 × 108 × 10 атомов и 972

фуллерена. В этом случае была рассмотрена только

Рис. 2. Расчет взаимодействия в приближении однородно-

кристаллическая пленка Al(111), а в расчетах, кото-

го вещества со средней плотностью. а — разбиение под-

рые проводились при T = 300 K, параметры потен-

ложки сферами радиуса L для интегрирования (3); б —

циала взаимодействия не занижались. В ходе вычис-

разбиение сферы параллельными окружностями для рас-

лений были определены скорости центров масс мо-

чета взаимодействия между атомом и однородной сферой;

в — разбиение подложки сферами для расчета взаимодей-

лекул в различные моменты времени. Полученные

ствия между сферой и подложкой с вакансией

результаты использовались для расчета спектраль-

ной плотности колебаний S(ω) по формуле [40]

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

∫∞

⁽∂V_f )

β =

(7)

S =2

〈v_z (t) · v_z (0)〉 cos(ωt) dt.

(4)

V_f

∂T

p=0

Индекс f в выражении (7) и ниже означает, что при

Также было рассчитано поверхностное натяже-

вычислении объема V_f использовано определенное

ние на границе Al-C₆₀ при температуре 1000 K для

ранее значение радиуса контактной поверхности R_f .

расплава алюминия и погруженных в него фул-

При погружении фуллерена в расплав под дейст-

леренов. Рассматривалась система из

61500

ато-

вием давления Лапласа p изменяется объем молеку-

мов: 54000 атомов алюминия и 7500 атомов угле-

лы на величину ΔV_f . По этим значениям, известным

рода (125 фуллеренов). Это соотношение компонен-

из МД, был рассчитан объемный модуль фуллерена:

тов примерно соответствует концентрации фуллере-

(

)

на 5.8 вес. %. В начальный момент времени атомы

∂p

K_T = -V_f

≈ -V_f

(8)

алюминия располагались в узлах ГЦК-решетки, а

∂V_fT

ΔV_f

фуллерены — в узлах простой кубической решетки.

Затем система стабилизировалась при 1000 K и ну-

Необходимо подчеркнуть некоторую неопреде-

левом давлении в течение 100 пс с использованием

ленность при вычислении объема, плотности, опре-

алгоритма, предложенного в работе [41]. После ста-

делении механических свойств и других макроско-

билизации термостат и баростат выключались и все

пических характеристик фуллерена, обусловленную

последующие расчеты проводились в микроканони-

малостью молекулы. Поскольку флуктуации мик-

ческом ансамбле. Временной интервал моделирова-

роскопических величин для малых молекул зна-

ния составлял 10 нс. Анализ проводился только для

чительны, параметры должны быть усреднены по

данных, полученных в микроканоническом ансамб-

ансамблю. В композиционных материалах такое

ле (NVE).

усреднение легко выполнить из-за большого количе-

В ходе моделирования вычислялись температу-

ства частиц в реальных образцах. Кроме того, зна-

ра, потенциальная энергия и локальное давление,

чения связанных с объемом или плотностью термо-

приходящееся на атом углерода. Для расчета ло-

динамических параметров β, γ или K_T чувствитель-

кального давления был использован алгоритм из ра-

ны к способу вычисления радиуса частицы. Расчеты

боты [42] и следующее выражение:

можно проводить с учетом или без учета конечного

объема атома углерода. Без указания метода расче-

n_V

Tr(mv_a · v_b + W_ab),

(5)

та радиуса фуллерена эти значения не имеют смыс-

ла.

где m — масса атома, v — его скорость, индексы a и

Влияние на результат выбора способа вычисле-

b соответствуют компонентам векторов, W — вири-

ния объемного модуля углеродных наночастиц, в

ал, Tr — след матрицы. Для расчета объемной плот-

том числе фуллеренов, подробно обсуждается в ста-

ности атомов углерода n_V предполагалось, что гра-

тье [43]. Показано, что применение модели, осно-

ница фуллерена совпадает со сферической поверх-

ванной на использовании для вычисления объемно-

ностью контакта Al/C₆₀, радиус которой задавал-

го модуля жесткости химической связи C-C, устра-

ся приближенно как R_f ≈ R + 2-5/6σ. Значение R

няет неопределенность, связанную с выбором объе-

вычислялось по данным МД как среднее расстоя-

ма наночастицы. Сжатие наночастиц при этом осу-

ние между атомами углерода и центром масс фул-

ществлялось методом перенормировки координат.

лерена. Результаты усреднялись по ансамблю и по

Заметим, что при условии сжатия частицы под дей-

времени в интервале 1 нс. Поверхностное натяжение

ствием внешней среды, которое обычно имеет ме-

рассчитывалось по давлению Лапласа:

сто в экспериментах, указанная проблема возника-

pR_f

ет вновь и проявляется в неопределенности упру-

γ =

(6)

гих свойств сжимающего вещества. Данный эффект

Далее в работе изучалась температурная зависи-

обусловлен применением континуальных моделей к

мость радиуса фуллерена R(T ) для свободных мо-

микроскопическим объектам и может быть нагляд-

лекул при нулевом давлении. При этом использова-

но проиллюстрирован при сравнении объемных мо-

лось усреднение по ансамблю из 1000 фуллеренов

дулей ГЦК-фуллерита [44], где межмолекулярное

и времени на 10000 шагов (2 пс). Рассматривался

взаимодействие фуллеренов осуществляется по дис-

температурный диапазон от 100 до 2000 К. Затем

персионному механизму и величина объемного мо-

вычислялся коэффициент температурного расшире-

дуля составляет 10.8 ГПа, и сверхтвердого фулле-

ния (КТР):

рита [45, 46], где фуллерены связаны ковалентно и

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

расчетные значения объемного модуля находятся в

Фуллерен также рассматривался как тонкая одно-

пределах от 236 ГПа до 304 ГПа. В первом слу-

родная сфера со средней плотностью поверхностно-

чае упругое сжатие вещества сопровождается пре-

го распределения атомов

имущественным изменением расстояний между со-

седними молекулами, тогда как во втором случае

n₂ =

≈ 0.3761Å-2.

4πR²

существенно также изменение внутримолекулярных

межатомных расстояний. Мы считаем возможным

Радиус сферы был рассчитан как среднее расстоя-

исследование упругих свойств фуллеренов в ука-

ние между атомами углерода и центром масс фул-

занных системах при рассмотрении средних по ан-

лерена после релаксации атомных позиций. Длины

самблю деформаций «отдельных» молекул под дей-

связи между атомами, расположенными в верши-

ствием внешней среды. При анализе зависимости от

нах пятиугольных структурных элементов фуллере-

давления радиусов молекул расчет объемных моду-

на, составляли 1.41Å, а между атомами в вершинах

лей фуллеренов с различными сжимающими веще-

шестиугольников — 1.45Å. Этот результат хорошо

ствами дал бы существенно меньшие различия, чем

согласуется с экспериментальными данными 1.40Å

полученные для кристаллов. Однако неоднозначно

и 1.46Å (см. [47] и цитированную там литературу).

определенными окажутся плотность и упругие мо-

Принятые допущения не позволяют учесть за-

дули сжимающего вещества, которое в этом слу-

висимость энергии связи от положения фуллерена

чае будет занимать дополнительный объем, прихо-

на поверхности алюминия, которая неявно пред-

дящийся на пространство между отдельными фул-

полагается слабой. Справедливость данной гипо-

леренами. Эту неопределенность можно устранить,

тезы неочевидна, но следующие эксперименталь-

включив в модель зависимость объемов сжимаемо-

но известные факты свидетельствуют в ее пользу.

го и сжимающего веществ от потенциала взаимодей-

Во-первых, согласно работам [9, 11] фуллерены ха-

ствия между атомами, относящимися к различным

рактеризуются высокой подвижностью на алюми-

молекулам.

ниевой подложке даже при комнатной температуре.

В настоящей работе предлагается приближен-

Во-вторых, температурная зависимость плотности

но учесть указанную неопределенность, рассчитав

покрытия n(T ) из экспериментов по ТПД [13] хо-

радиус R_f тонкой сферической границы между

рошо аппроксимируется аналитической функцией с

алюминием и фуллереном, причем R_f определяет-

единственной точкой перегиба. Поэтому мы предпо-

ся параметрами потенциала взаимодействия атомов

лагаем, что изменения энергии связи при смещении

Al-C. Считая величину R_f радиусом включения и

фуллеренов в плоскости не играют существенной ро-

используя усреднение по ансамблю, можно полу-

ли в экспериментах.

чить хорошо определенные для системы Al-C₆₀ зна-

Использование приближения однородной среды

чения термодинамических величин как для включе-

вместе с потенциалом межатомного взаимодействия

ния, так и для матрицы. Отметим, что в данной ра-

(1) позволяет получить формулу (2) для взаимодей-

боте значение R_f было вычислено в нулевом прибли-

ствия между атомом и подложкой. Аналогично, ин-

жении, учитывающем зависимость R_f от положения

тегрируя по поверхности сферы (см. рис. 2б), мы

минимума, но без учета жесткости потенциала (1).

получаем следующее выражение для энергии взаи-

модействия атома и фуллерена:

[

]

3. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ И

ϵ = 8πεσ²

(ζ-10 - ζ⁺¹⁰) -

(ζ-4 - ζ⁺⁴) ,

(9)

ОБСУЖДЕНИЕ

n₂ H

3.1. Взаимодействие фуллеренов с

где для сокращения записи использовано обозначе-

подложкой

ние

ζ^±m = [σ/(H ± R)]^m,

Для теоретического анализа численных резуль-

татов моделирования взаимодействия фуллеренов

H — расстояние между атомом и центром сферы,

с алюминиевой подложкой использовались следую-

R — ее радиус, m — показатель степени.

щие упрощающие допущения. Считалось, что атомы

Предположим, что вакансия имеет форму сфе-

алюминия распределены в подложке-полупростран-

рического сегмента, образованного пересечением

стве равномерно с заданной средней плотностью n₁.

алюминиевой подложки и сферы, центр которой

Значение n₁ ≈ 0.0590Å-3 было рассчитано после

совпадает с центром масс фуллерена (см. рис. 2в).

оптимизации элементарной ячейки в рамках DFT.

Параметр L₀ определим из равенства объемов

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

удаленного сегмента и атома алюминия. При этом

H = 5.50Å. При наличии вакансии на поверхности

будем считать, что расположение фуллерена над

параметры существенно отличаются: W₀ = -2.04 эВ

подложкой соответствует минимуму потенциальной

и H = 4.99Å. Этот эффект соответствует описан-

энергии H = H₀:

ному в работе [14]. В качестве одной из причин в

[

]

статье [14] был отмечен сильный стерический эф-

n-11 = π

(L³⁰ - H³⁰) - H₀(L²⁰ - H²⁰)

(10)

фект, наблюдаемый при релаксации атомных пози-

ций. В ходе оптимизации ближайший к фуллерену

Интегрируя (9) по объему алюминия, получа-

атом алюминия смещался вглубь относительно по-

ем зависимость потенциальной энергии от рассто-

верхности подложки. Смещение атома было вызва-

яния H:

но чрезмерным сближением атомов алюминия и уг-

лерода и сопровождалось увеличением расстояния

W (H, L₀) =

π²εRn₁n₂σ³ ×

H и энергии E_Al-C. Из геометрических соображений

[

]

очевидно, что удаление указанного атома алюми-

(I-9 - I⁺⁹) - (I-3 - I⁺³) ,

(11)

ния приводит к увеличению прочности связи меж-

ду фуллереном и подложкой и смещению молекулы

где введено обозначение

к поверхности. Кроме того, эффект усиления свя-

I^±m = [mL₀ - (m - 1)H ± R]ζ^±m,

зи может быть вызван увеличением химической ак-

тивности атомов алюминия за счет уменьшения их

m принимает значения 3 или 9. В отсутствие вакан-

координационного числа.

сии L₀ = H и выражение (11) сводится к виду

Использование эмпирической модели межатом-

ного взаимодействия (1) позволяет разделить воз-

W (H) =

π²εRn₁n₂σ⁴ ×

можные причины этого эффекта. С одной сторо-

[

]

ны, поскольку параметры потенциала не зависят от

(ζ-8 - ζ⁺⁸) - (ζ-2 - ζ⁺²)

(12)

относительного положения атомов алюминия, мо-

дель нечувствительна к возможному изменению хи-

Вычисленные в рамках DFT значения H = H₀ и

мической активности атомов при возникновении ва-

W (H₀) = W₀ соответствуют минимуму потенциаль-

кансии. С другой стороны, потенциал (1) учиты-

ной энергии. Применив условие минимума потенци-

вает отталкивание атомов на малых расстояниях

альной энергии к выражению (12) и приравнивая

и позволяет воспроизвести особенности взаимодей-

полученные при этом значения H₀ и W₀ величи-

ствия, связанные со стерическим эффектом. Кроме

нам, предварительно вычисленным ab initio, запи-

того, выражение (11) дает возможность учета из-

шем систему уравнений для определения парамет-

менения энергии взаимодействия и смещения фул-

ров потенциала (1):

лерена. Выражение (11) было получено в предпо-

]1/6

ложении о специфической форме вакансии с цент-

^[15 (H₀ + R)³ - (H₀ - R)

σ = (H²⁰-R²)

(13)

ром кривизны, совпадающим с центром фуллерена

(H₀ + R)⁹ - (H₀ - R)⁹

(см. рис. 2в). Это допущение существенно упрощает

процедуру интегрирования. Можно предположить,

[

]-1

3W₀

что использование более реалистичной модели поз-

ε=

(ζ-8 - ζ⁺⁸)-(ζ-2 - ζ⁺²)

(14)

2π²n₁n₂Rσ⁴

волит точнее воспроизвести результаты расчетов ab

initio. Из-за специфической формы вакансии, соот-

Энергия связи фуллерена с подложкой W₀ рас-

ветствующей форме поверхности фуллерена, анали-

считывалась с использованием данных DFT по сле-

тическая оценка с помощью формулы (11) может

дующей формуле:

завышать абсолютное значение энергии связи и за-

нижать равновесное расстояние H₀. Тем не менее мы

W₀ = E_Al-C - (EAl + EC),

(15)

полагаем, что при условии дополнения результата-

где E_C и E_Al — потенциальные энергии свободной

ми более строгих численных расчетов предложенная

молекулы и алюминиевой подложки соответственно,

упрощенная модель может быть полезна для анали-

за обсуждаемых эффектов.

а E_Al-C — энергия фуллерена, взаимодействующего

с поверхностью подложки.

Энергия взаимодействия между фуллеренами и

Расчеты в рамках DFT для фуллерена, распо-

подложкой была рассчитана при нулевой темпера-

ложенного на идеальной поверхности Al(111), дают

туре с использованием как модели ab initio, так и

значение энергии связи W₀ = -0.98 эВ и расстояние

эмпирического потенциала вида (1) с параметра-

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

Таблица. Взаимодействие фуллерена с поверхностью Al(111) при нулевой температуре. Расчет с использованием

DFT и потенциала (1)

W₀, эВ

H₀,Å

C-Al,Å

Wvac0, эВ

Hvac0,Å

C-Al^vac,Å

DFT (эта работа)

-0.98

5.50

2.24

-2.04

4.99

2.19

DFT [14]

-1.37

2.24

-2.34

2.20

Эмпирический

-1.86

5.20

2.31

-1.96

4.88

2.25

потенциал (1)

ми, определенными согласно (13) и (14). Результа-

Выражения (11) и (12) дают основание предпо-

ты в сравнении с данными работы [14] представ-

лагать, что изменение параметра ε с постоянной

лены в таблице. Для удобства сравнения с рабо-

σ не должно влиять на значение H₀, а W₀ ли-

той [14] в таблице добавлены усредненные значе-

нейно зависит от ε. Однако оптимизация с помо-

ния длин связей C-Al. Оптимизация положения ато-

щью программы LAMMPS свидетельствует об об-

мов и расчет энергии с использованием потенци-

ратном. Причина расхождения в том, что при вы-

ала (1) в программе LAMMPS позволяет учиты-

воде уравнений (11) и (12) подразумевалась идеаль-

вать фактическое расположение частиц без исполь-

ная жесткость подложки и сферы, в численных же

зования каких-либо предположений об однородно-

расчетах использовались реалистичные потенциалы

сти вещества. Сравнение результатов эмпирических

для атомов углерода и алюминия [35,37], учитыва-

и ab initio расчетов позволяет нам оценить приме-

ющие смещение атомов и деформацию поверхнос-

нимость модели однородного вещества и выражений

тей. Подставляя в выражения (13) и (14) значе-

(13), (14) для параметризации потенциала (1).

ния Wvac0 = -2.04 эВ и H = 5.50Å, рассчитанные

Результаты вычислений ab initio настоящей ра-

в настоящей работе, получаем ε = 5.27 · 10-2 эВ,

боты находятся в разумном соответствии с данны-

σ = 2.70Å. Оптимизация в программе LAMMPS с

ми [14]. Некоторое расхождение в значениях энергий

использованием указанных выше параметров дает

H₀ = 5.3Å для случая без вакансии и Hvac0 = 4.9Å

связи W₀ и Wvac0 может быть связано с различиями

в используемых расчетных моделях.

для случая с вакансией. Такое изменение положе-

ния фуллерена хорошо согласуется с результатами

Расчет параметров потенциала по формулам (13)

расчетов в рамках DFT. Значение W₀ при расчете

и (14) для подложки без вакансии дает значения ε =

с использованием потенциала (1) изменяется намно-

= 2.67 · 10-2 эВ, σ = 2.70Å. Подставляя эти пара-

го меньше, чем значение, вычисленное ab initio. Без

метры в выражение (11), получаем для подложки с

вакансии энергия связи, рассчитанная в программе

вакансией Wvac0 = -1.46 эВ, Hvac0 = 4.77Å. Следова-

LAMMPS, составила W₀ = -1.86 эВ, а с ваканси-

тельно, предложенная аналитическая модель лишь

ей Wvac0 = -1.96 эВ. Следовательно, при исполь-

частично воспроизводит эффекты, связанные с из-

зовании потенциала (1) с предложенными парамет-

менением энергии W₀ и расстояния H₀ из-за обра-

рами для расчета взаимодействия C₆₀ с идеальной

зования вакансии на поверхности подложки. Для

подложкой Al(111) энергия связи оказывается силь-

расстояния H₀ аналитическая модель дает несколь-

но завышенной. Мы полагаем, что это несуществен-

ко заниженное значение по сравнению с DFT, что

но для исследования композиционных материалов

может быть связано с использованием упрощенной

Al-C₆₀, поскольку технологии производства объем-

формы вакансии. Энергия связи по абсолютной ве-

ных композитов предполагают высокоэнергетиче-

личине сильно занижена (1.46 эВ против 2.04 эВ).

скую обработку вещества, активирующую процес-

Расчет с использованием программы LAMMPS с

сы, связанные с реконструкцией поверхности алю-

указанным эмпирическим силовым полем не пред-

миния.

сказывает существенного изменения W₀ и H₀ из-за

вакансии. Как с вакансией, так и без нее оптимиза-

Подстановка в формулы (13) и (14) при вычис-

ция дает значения W₀ ≈ -0.9 эВ и H₀ ≈ 5.4Å. Сле-

лении параметров Wvac0 = -2.04 эВ и Hvac0 = 4.99Å

довательно, модель дает правильные значения для

дает значения ε = 0.175 эВ и σ = 1.99Å. Исполь-

идеальной поверхности Al(111), но не предсказыва-

зование таких параметров при оптимизации пози-

ет изменений, связанных с появлением вакансии.

ций атомов в программе LAMMPS приводит к зна-

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

чительному искажению структуры подложки и к

Для разработки аналитической модели были ис-

неправдоподобным значениям W₀ и H₀.

пользованы следующие упрощающие приближения.

Проведенный анализ позволяет сделать вывод,

Фуллерены рассматривались как точечные части-

что наилучшее соответствие с результатами расче-

цы в пренебрежении внутренними степенями сво-

тов DFT дает потенциал (1) с параметрами ε =

боды и вращательной энергией. Движение частиц

= 5.27 · 10-2 эВ, σ = 2.70Å. Поэтому далее в на-

в плоскости предполагалось свободным и неогра-

стоящем исследовании были использованы эти зна-

ниченным. Это предположение следует из экспери-

чения. Сравнение теоретических расчетов с экспе-

ментов [9,11], где высокая подвижность фуллеренов

риментальными данными [9,13] также указывает на

в плоскости подложки наблюдалась даже при ком-

предпочтительное использование выбранных значе-

натной температуре. Кинетика десорбции исследо-

ний параметров (см. разд. 3.2). Данные, представ-

валась с помощью теории переходного состояния.

ленные в таблице, также были получены с использо-

Свободная энергия активации вычислялась в гармо-

ванием этих значений параметров. Из таблицы вид-

ническом приближении и принималась равной раз-

но, что результаты, полученные с использованием

ности свободных энергий продуктов реакции и ре-

потенциала (1) для поверхности алюминия с вакан-

агентов. В пренебрежении адсорбцией предполага-

сией, лучше согласуются с расчетами ab initio. От-

лось, что удаленные с поверхности фуллерены боль-

метим, что потенциал (1) с заданными параметрами

ше не могут взаимодействовать с подложкой.

существенно завышает энергию связи для фуллере-

Принятые упрощения позволяют рассчитать ско-

нов и идеальной подложки Al(111). Основной при-

рость процесса десорбции, используя следующее

чиной различий в энергиях связи W₀ и Wvac0 являет-

уравнение [16]:

ся зависимость химической активности атомов алю-

(

)

nω

ΔF

миния и углерода от их локального окружения, ко-

k=-

exp

(17)

2π

k_BT

торая не учитывается потенциалом (1). Изменение

положения фуллерена относительно подложки свя-

где k = dn/dt — скорость десорбции, n — плотность

зано с конечностью объемов атомов и при исполь-

покрытия (количество молекул на единице площа-

зовании потенциала (1) воспроизводится с высокой

ди), в начальный момент времени равная n₀, ΔF —

точностью.

свободная энергия активации, ω — частота колеба-

Заметим, что использование парного потенциа-

ний молекулы в направлении нормали к поверхнос-

ла (1) не позволяет учесть имеющее место перерас-

ти, k_B — постоянная Больцмана.

пределение заряда при формировании полярной ко-

Использование распределения Больцмана в вы-

валентной связи Al-C. Также не учитывается из-

ражении (17) накладывает определенные ограниче-

менение параметров взаимодействия для пар ато-

ния на условия задачи. Время релаксации десорб-

мов Al-Al и C-C, обусловленное перераспределени-

ции должно быть намного больше, чем время темпе-

ем электронной плотности при формировании свя-

ратурной релаксации системы, возмущенной актом

зи Al-C и выступающее причиной изменения ло-

элементарной реакции. В настоящей модели это тре-

кальной структуры фуллерена [14]. Использование

бование предполагается выполненным.

же существующих моделей, приближенно учитыва-

Решение уравнения (17) имеет следующий общий

ющих перенос заряда при образовании химической

вид:

связи и зависимость ее энергии от локального окру-

⎡

⎤

∫

(

)

жения атомов (например, [48,49]), ограничено ввиду

ΔF

n(t) = n₀ exp^⎣-

exp

dt^⎦ .

(18)

отсутствия параметризации с включенными в тре-

2π

k_BT

нировочный набор системами Al-C₆₀.

Используя гармоническое приближение, можно

аналитически вычислить конфигурационный инте-

3.2. Кинетика десорбции фуллеренов с

грал и получить выражение для свободной энергии

подложки при линейном изменении

активации:

температуры во времени

⁽2πk_BT⁾

Рассмотрим процесс десорбции фуллеренов с по-

ΔF = -W₀ +

k_BT ln

(19)

ασ²

верхности вещества при нагревании системы от тем-

пературы T₀ до T₁ с постоянной скоростью c:

Здесь жесткость α определяется значением второй

производной потенциальной энергии (12) и связана

T = T₀(1 + ct).

(16)

с частотой колебаний и массой осциллятора:

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

)1/2

^√α

⁽1 d²W

ω=

(20)

m dH²

Подставляя W (H) из (12) в (20), получим

[

]

α = 4π²n₁n₂Rσ²ε

(ζ-10 - ζ⁺¹⁰) - (ζ-4 - ζ⁺⁴)

(21)

В принятой постановке в подынтегральном вы-

ражении правой части (18) от времени зависит толь-

ко температура. В статьях [9, 11] отжиг проводился

при постоянной температуре. В этом случае инте-

грал под знаком экспоненты дает линейную функ-

цию времени, и эволюция плотности покрытия опре-

деляется выражением

n(t) = n₀e^-t/τ ,

(22)

где время релаксации τ — константа,

Рис. 3. Спектральная плотность колебаний центров масс

(

)-1/2

(

)

ασ

фуллеренов в направлении нормали к алюминиевой под-

τ = 2πω-1

exp

(23)

2πk_BT

k_BT

ложке

Рассмотрим случай линейного нагрева системы,

реализованный в работе [13]. Тогда после интегри-

которое, в свою очередь, совпадает с формулами

рования получим следующие выражения для плот-

(22) и (23) с точностью до малой поправки ко вре-

ности покрытия:

мени релаксации δτ ∝ 3k_BT₀/(2W₀).

С учетом данных из таблицы и рассчитанных

n≈n₀×

по ним значений параметров ε и σ были вычисле-

{

√

}

W₀ασ²

ны жесткость α и частота ω согласно (20) и (21):

× exp

[ψ(T ) - ψ(T₀)]

(24)

α = 9.11 эВ·Å-2, ω = 10.8 пс-1. Для проверки ана-

2πc

^-2πk2B

T20

литических расчетов и анализа применимости моде-

√

ли был проведен спектральный анализ траекторий

(

)

nω ασ

W₀

центров масс фуллеренов. Расчетная спектральная

exp

(25)

2π

2πk_BT

k_BT

плотность колебаний представлена на рис. 3.

Максимум функции S(ω) соответствует частоте

В выражении (24) мы ограничились первыми дву-

ω = 9.1 пс-1, что в пределах 20% совпадает с ана-

мя членами асимптотического разложения дополни-

литическим значением. Принимая во внимание зна-

тельной функции ошибок erfc(x) для больших аргу-

чимость принятых предположений, мы считаем та-

ментов x = |W₀|/(k_B T ) ≫ 1 и использовали следу-

кое совпадение результатов допустимым. Использо-

ющее обозначение:

вание слегка завышенной частоты колебаний в ана-

(

)-3/2

(

)

литических расчетах может привести к переоценке

W₀

ψ(T ) =

exp

(26)

скорости десорбции.

k_BT

Расчет времени релаксации десорбции по фор-

Используя разложение в ряд Тейлора, можно по-

муле (23) при температуре 730 К дает τ ≈ 5.3 с.

казать, что в пределе c → 0 уравнение (24) сводится

Используя значение из работы [14] W₀ = 2.34 эВ, по-

к выражению

лучим τ ≈ 544 с. Учитывая условия эксперимента в

⎡

√

[9,11], заметим, что полученные времена релаксации

(

)3/2

ασ

W₀

вполне реалистичны. Можно предположить, что ре-

n ≈ n0 exp⎣ωt

2π

2πk_BT₀

k_BT₀

зультаты работы [14] несколько завышают абсолют-

⎤

ную величину энергии связи, а данные настоящей

(

)

работы — занижают. Использование значения энер-

3k_BT₀

× ψ(T₀)

⎦,

(27)

гии связи W₀ = -0.98 эВ, полученного для идеаль-

2W₀

ной подложки Al(111), ведет к существенному зани-

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

Рис. 4. Температурные зависимости плотности покрытия фуллеренами алюминиевых подложек. Сплошные линии — ана-

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

для жидкой пленки алюминия и W0 = -1.41 эВ; б — десорбция с кристаллической подложки Al(111) при W0 = -0.51 эВ

жению времени релаксации τ ≈ 5.3 · 10-7 с. Анализ

чение параметра σ не зависит от W₀, поэтому счи-

показывает, что при такой энергии связи десорбция

талось, что σ = 2.70Å.

фуллеренов с поверхности алюминия будет наблю-

Температурные зависимости плотности покры-

даться в результате отжига при температуре 350 K

тий, рассчитанные аналитически и методом МД в

со временем релаксации 8 с, что противоречит экс-

указанных условиях, представлены на рис. 4. Хоро-

периментальным данным [9, 11].

шее согласие численных и аналитических результа-

Расчет методом МД ограничен временами поряд-

тов свидетельствует об адекватности модели. Ана-

ка 1 нс, а время десорбции с параметрами, опреде-

литический расчет несколько занижает температу-

ленными выше, намного больше. Формально, умень-

ру десорбции по сравнению с численным. Частично

шив значение энергии связи W₀ и увеличив темпе-

этот эффект связан с рассмотренной ранее неточно-

ратуру, мы можем значительно ускорить десорбцию.

стью при вычислении частот колебаний. Кроме то-

Несмотря на то, что условия модели при этом не со-

го, отметим, что расхождение более выражено для

ответствуют экспериментальным, сравнение анали-

расплава, чем для кристаллической подложки. Уси-

тических результатов с данными прямого МД-моде-

ление взаимодействия фуллеренов с поверхностью

лирования необходимо для проверки как аналити-

жидкости может быть связано с капиллярными яв-

ческой модели, так и применимости потенциала (1).

лениями.

Моделирование методом молекулярной динами-

Моделирование десорбции методом МД проводи-

ки не предсказывает образования карбида Al₄C₃

лось для различных температур и параметров по-

при высоких температурах, которое наблюдалось в

тенциала взаимодействия. В первом случае темпе-

экспериментах с фуллеренами и углеродными нано-

ратура системы повышалась от 1000 K до 2200 K,

трубками [50, 51]. Это объясняется малыми време-

а значение энергии связи задавалось равным 70 %

нами МД-моделирования и активационным харак-

от рассчитанного в рамках DFT, W₀ = -1.43 эВ.

тером химической реакции.

Этот диапазон температур соответствует равнове-

сию алюминия в жидкой фазе. Во втором слу-

В условиях исследований методом ТПД скорос-

чае рассматривалась кристаллическая подложка

ти нагрева могут варьироваться в широком диапа-

Al(111), а нагрев проводился в температурном диа-

зоне значений [16]. Для анализа кинетики десорбции

пазоне от 400 K до 800 K. При этом энергия свя-

в настоящей работе рассмотрены скорости нагрева

зи была уменьшена в четыре раза по сравнению с

cT₀ = 0.1, 3, 10 К/с (обозначения переменных со-

расчетной, W₀ = -0.51 эВ. В обоих случаях время

ответствуют (16)). Теоретические зависимости n(T)

нагрева задавалось равным 1 нс. Согласно (13), зна-

сравнивались с экспериментальными данными [13],

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

Рис. 5. Температурные зависимости плотностей покрытия поверхности (a) и скоростей десорбции (б), рассчитанные по

уравнениям (24) и (25). Сплошные линии соответствуют скорости нагрева 0.1 К/с, штрихпунктирные линии — 3 К/с,

штриховые линии — 10 К/с. Экспериментальные данные из [13], полученные при скорости нагрева 3 K/с, обозначе-

ны крестиками. Пунктирной линией показаны результаты аналитических расчетов, выполненных в настоящей работе по

уравнению (24) со значением W0 = -2.34 эВ из [14]

где скорость нагрева составляла 3 К/с. Результаты

3.3. Капиллярные явления на границах

представлены на рис. 5a. Кроме того, предложенная

между расплавом алюминия и

модель была использована для расчета кинетики де-

погруженными в него фуллеренами C₆₀

сорбции при скорости нагрева 3 K/с, соответствую-

щей эксперименту [13], для значения W₀ = -2.34 эВ,

Анализ состояния расплава алюминия с погру-

взятого из [14]. Полученный результат сравнивался

женными в него фуллеренами указывает на нерав-

с данными настоящей работы и экспериментом [13]

новесность системы. В этом случае можно выделить

(см. рис. 5a). Расчетные температурные зависимо-

несколько процессов с разными характерными вре-

сти скорости десорбции при различных скоростях

менами. Релаксация функции распределения ато-

нагрева показаны на рис. 5б.

мов по скоростям в данном случае является наибо-

Из рис. 5 видно, что при использовании резуль-

лее быстрым процессом, характерное время которо-

татов расчетов ab initio из настоящей статьи предло-

го составляет τ_at ∼ 1-10 пс, т. е. в 10-100 раз больше

женная аналитическая модель воспроизводит экспе-

периода колебаний атомов. Установление равновес-

риментальные данные [13] с высокой точностью. Ис-

ного поступательного движения фуллеренов зани-

пользование значения энергии связи из работы [14]

мает немного больше времени, τC60 ∼ 5τ_at. Прибли-

приводит к замедлению десорбции и ухудшению со-

женно время этого процесса можно оценить, учиты-

гласия между теорией и экспериментом. С ускоре-

вая, что при лобовом упругом столкновении частиц

нием нагрева максимум скорости десорбции смеща-

с массами m₁ и m₂ (m₁ < m₂) величина перераспре-

√

ется в область более высоких температур. Максиму-

деленной энергии ΔE ∝

m₁/m₂. Релаксация про-

мы скорости десорбции при скоростях нагрева 0.1,

странственного распределения фуллеренов опреде-

3 и 10 К/с соответствуют температурам 685, 757 и

ляется коэффициентом диффузии D наночастиц в

785 К. Полученный диапазон от 685 K до 785 K пол-

расплаве. Значение D ∼ 10-5 см²/с может быть при-

ностью согласуется с данными [9,11], согласно кото-

нято в качестве разумной оценки. Время релаксации

рым температура десорбции составляет 730 K. Ре-

процесса диффузии τ_D = 1/(4Dns/30) ∼ 10-10 с, где

зультаты моделирования указывают на единствен-

ns0

≈ 1.15 · 10-4 Å-3 — концентрация фуллеренов

ный пик функции k(T ), что является следствием ис-

в единице объема. Кроме того, из теории коллоид-

пользованного предположения о постоянстве энер-

ных растворов известно, что при слабой смачивае-

гии связи W₀, независимо от расположения фулле-

мости окружающей жидкостью частицы склонны к

ренов на подложке.

коагуляции под действием капиллярных сил [52-54].

ЖЭТФ, вып. 1

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

Процесс кластеризации фуллеренов в расплаве алю-

миния наблюдался и в наших расчетах. Для пред-

варительной оценки характерного времени процес-

са воспользуемся допущениями, принятыми в тео-

рии Смолуховского [53, 54]. Предположим, что ан-

самбль фуллеренов можно представить как идеаль-

ный газ броуновских частиц. Мы также предпола-

гаем, что элементарный акт коагуляции происходит

всегда при сближении двух фуллеренов на расстоя-

ние, меньшее или равное 2R_f ≈ 10.4Å. Для описа-

ния процесса используем кинетическое уравнение:

dn_s

= -k_Cn2s,

(28)

где n_s

— концентрация свободных фуллеренов,

которая рассчитывалась непосредственно в ходе

МД-моделирования. В начальный момент времени

Рис. 6. Зависимость обратной концентрации свободных

n_s = ns0. Константа скорости коагуляции k_C в тео-

фуллеренов от времени: кружками обозначены результа-

рии Смолуховского определяется выражением

ты расчетов МД, сплошная линия — аппроксимация вы-

ражением (30), где kC ≈ 4.84·10-13 см³· с-1 — константа

k_C = 8πR_fD.

(29)

скорости коагуляции

Из уравнения (28) следует линейная зависимость

обратной концентрации свободных фуллеренов от

времени:

(30). Связанными при этом считались фуллерены,

ns0

= 1 + ns0k_Ct.

(30)

расстояние между центрами масс которых не превы-

n_s

шает 2(R+R_c) ≈ 9.5Å (R_c = 2.1Å — радиус обрезки

В качестве характерного времени коагуляции

потенциала Терсоффа [35]).

обычно используется величина τ1/2 — время, за ко-

На рис. 6 представлена зависимость обратной

торое концентрация свободных частиц уменьшается

концентрации свободных молекул от времени.

в два раза. Используя выражения (29) и (30), можно

Характерное время процесса коагуляции τ1/2

оценить время коагуляции как

= (k_C ns0)-1

≈ 18.0 нс, что примерно в 25 раз

больше значения, полученного с использованием

τ1/2 = (8πns0 R_f D)-1 ≈ 7 · 10-10 с.

теории Смолуховского. На временных интервалах

t₀ ≪ τ1/2 процессом коагуляции можно пренебречь.

Выполненные оценки масштабов времен релак-

В частности, за время t₀

= 1 нс после начала

сационных процессов предполагают разумный вы-

вычислений около 5 % фуллеренов образовали свя-

бор времени МД-моделирования 10 нс. Согласно

занные кластеры. Следовательно, для временного

оценкам, использующим предположение об отсут-

интервала t ≤ 1 нс можно рассчитать коэффициент

ствии взаимодействия между фуллеренами в рас-

диффузии свободных фуллеренов в расплаве. На

плаве, времена диффузии и коагуляции являются

рис.

показана зависимость среднеквадратично-

величинами одного порядка малости, т. е. невозмож-

го смещения (MSD) фуллеренов от времени при

но разделить эти процессы на «быстрые» и «медлен-

t ≤ 1 нс.

ные». Однако результаты моделирования методом

МД свидетельствуют о различных скоростях диф-

В простых жидкостях MSD имеет два характер-

фузии и коагуляции и позволяют сделать вывод о

ных режима, для разделения которых удобно ис-

сильном взаимодействии между фуллеренами.

пользовать график в логарифмическом масштабе.

В настоящей работе для более точного расчета

В баллистическом режиме зависимость имеет вид

скорости коагуляции был выполнен анализ траек-

〈Δr²〉 = v²t², а в диффузионном режиме — 〈Δr²〉 =

торий центров масс фуллеренов, вычисленных мето-

= 6Dt. Из рис. 7 видно, что диффузионный режим

дом МД. Константа скорости коагуляции была опре-

соответствует времени t > 0.2 нс. Линейная аппрок-

делена путем аппроксимации временной зависимо-

симация численных результатов выполнялась для

сти концентрации свободных фуллеренов формулой

этого временного интервала. Полученный в резуль-

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

Рис. 7. Зависимости MSD от времени в линейном (а) и логарифмическом (б) масштабах. Кружками обозначены резуль-

таты вычислений МД, сплошные линии — линейная аппроксимация

тате вычислений методом МД коэффициент диффу-

зии D = 3.64 · 10-5 см²/с по порядку величины со-

ответствует использованному ранее в оценках зна-

чению D ∼ 10-5 см²/с. Результаты расчетов поз-

воляют уточнить время диффузионной релаксации,

t_D = 0.29 нс. Пересчет константы скорости коагу-

ляции по формуле Смолуховского (29) дает k_C ≈

≈ 4.77 · 10-11 см³· с-1, что соответствует характер-

ному времени τ1/2 ≈ 0.18 нс.

Сравнение времен релаксации подтверждает воз-

можность разделения неравновесных процессов коа-

гуляции и диффузии на «медленные» и «быстрые».

Подчеркнем, что оценка константы скорости коагу-

ляции согласно теории Смолуховского, предполага-

ющей отсутствие взаимодействия между фуллере-

нами, дает значение примерно в 100 раз превышаю-

Рис. 8. RDF центров масс фуллеренов, погруженных в

щее результат, полученный из МД-моделирования.

расплав алюминия. Сплошная линия соответствует RDF,

Поэтому можно сделать вывод, что взаимодействие

усредненной на временном интервале t ∈ [0, 1] нс, штрихо-

между фуллеренами отлично от нуля на расстояни-

вая линия — t ∈ [5, 6] нс, пунктирная линия — t ∈ [9, 10] нс.

ях, превышающих 2(R + R_c), а процесс коагуляции

Вставка иллюстрирует расположение фуллеренов в рас-

носит активационный характер.

плаве, соответствующее первым двум пикам

На рис. 8 показаны функции радиального рас-

пределения (RDF) центров масс молекул C₆₀. Функ-

ции усреднялись по разным временным интервалам:

моделирования для взаимодействия между атомами

t ∈ [0,1] нс, t ∈ [5,6] нс и t ∈ [9,10] нс.

углерода был использован только короткодейству-

Эволюция RDF связана с ростом узкого пика в

ющий потенциал: радиус обрезки потенциала Тер-

окрестности точки rn1/3 = 0.45 (9.24Å), который ха-

соффа [35] равен 2.1Å. В случае, если расстояние

рактеризует количество связанных молекул. Функ-

между центрами масс больше 9.5Å, потенциал Тер-

ция g(r) имеет ярко выраженные максимумы в точ-

соффа не вносит вклада в межмолекулярное взаи-

ках 0.71, 0.82 и 0.94 (14.6, 16.8 и 19.3Å), что указы-

модействие. Можно предположить, что взаимодей-

вает на неидеальность ансамбля фуллеренов. В ходе

ствие фуллеренов на больших расстояниях связано с

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

упругой деформацией расплава алюминия под дей-

ствием капиллярных сил.

Подчеркнем, что в настоящей работе не учитыва-

лось имеющее место перераспределение заряда при

образовании полярных ковалентных связей Al-C, и

фуллерены считались электрически нейтральными.

Мы предполагаем, что учет кулоновского отталки-

вания между фуллеренами приведет к замедлению

процесса коагуляции, поэтому рассчитанное в насто-

ящей работе значение константы скорости можно

считать оценкой сверху. Исследование влияния ку-

лоновских поправок к потенциалу межмолекуляр-

ного взаимодействия на скорость коагуляции фул-

леренов в расплаве выходит за рамки настоящей ра-

боты.

Изменение параметров состояния в результа-

те релаксации системы при образовании класте-

Рис.

9. Температурная зависимость объема фуллере-

ров связанных фуллеренов несущественно. Разница

на. Кружками обозначены результаты МД-моделирования,

между значениями температуры и давления, рас-

сплошной линией — аппроксимация V_f /V0 = βT + 1.

считанными и усредненными за первую и десятую

V0 = 589.79Å3 — объем фуллерена при 0 K, рассчитанный

наносекунды, составляет 3 K и 300 бар, а сред-

путем экстраполяции результатов вычислений, β — КТР

неквадратические отклонения составляют 0.85 K и

80 бар. Таким образом, изменение параметров со-

стояния во время релаксации по порядку величины

Рассчитанный по формуле (7) с использованием

совпадает со статистической ошибкой.

линейной аппроксимации результатов МД КТР со-

Расчет давления Лапласа, действующего на фул-

ставил β = 1.15 · 10-5 K-1. Представленные резуль-

лерены со стороны расплава алюминия, дает значе-

таты отличаются от данных [56], где авторы получи-

ние p = 1.49 ГПа. Эта величина медленно изменяет-

ли нелинейную зависимость V (T ) в диапазоне тем-

ся при коагуляции фуллеренов: изменение за 10 нс

ператур от 0 до 400 K, а КТР принимал значения

в ходе моделирования не превысило 0.5 %. Положи-

от -1 · 10-5 K-1 до 2.0 · 10-5 K-1. Отметим, что ав-

тельное значение p указывает на слабую смачива-

торами [56] получена сильная зависимость β(T) при

емость фуллеренов расплавом. Поверхностное на-

низких температурах. При повышении температу-

тяжение, рассчитанное по формуле (6), составляет

ры выше 300 К КТР медленно возрастает, асимп-

γ = 0.39 Дж/м². Согласно [55], коэффициент по-

тотически приближаясь к некоторому постоянному

верхностного натяжения алюминия, рассчитанный

значению β ≈ 2.0 · 10-5 K-1. Значение КТР, рас-

с тем же потенциалом при 1000 К, составляет при-

считанное в настоящей статье, заметно меньше, чем

мерно 0.71 Дж/м², т. е. больше на 0.32 Дж/м², а

рассчитанное в [56], что иллюстрирует уже обсуж-

для межфазной границы между расплавом алюми-

давшийся ранее факт сильной зависимости парамет-

ния и графитовой подложкой — 0.18 Дж/м², т. е.

ров β, γ и K от метода расчета радиуса фуллерена.

меньше на 0.21 Дж/м². Поскольку притяжение меж-

Напомним, что в настоящей работе использован ра-

ду атомами углерода и алюминия в случае фулле-

диус контактной поверхности R_f = R + 2-5/6σ, а в

ренов намного сильнее, чем для графита, несколь-

[56] авторы использовали среднее расстояние от ато-

ко неожиданно, что свободная энергия единицы по-

мов до центра масс фуллерена R. Пересчет КТР с

верхности Al/C₆₀ оказывается выше. По-видимому,

использованием R дает значение β = 1.63·10-5 K-1,

данный эффект связан со сферической геометрией

что неплохо согласуется с данными [56] при темпе-

фуллерена и меньшим количеством участвующих во

ратурах выше 250 К.

взаимодействии атомов.

Среди возможных причин расхождения резуль-

Для определения КТР фуллеренов анализиро-

татов при более низких температурах можно выде-

вался ансамбль из 1000 невзаимодействующих мо-

лить использование потенциала Терсоффа для опи-

лекул. Расчетная зависимость объема фуллерена от

сания взаимодействия атомов углерода. Согласно

температуры близка к линейной (см. рис. 9).

результатам [56], где расчет сил и энергии выпол-

100

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

нялся в рамках DFT, потенциал Терсоффа неудо-

ми C₆₀. С использованием DFT рассчитана энергия

влетворительно описывает температурную зависи-

связи фуллеренов с подложкой Al(111). Определены

мость объема углеродных наночастиц при низких

параметры потенциала Леннарда-Джонса для пар

температурах.

атомов Al-C. Предложена аналитическая модель де-

Под действием давления Лапласа p = 1.49 ГПа

сорбции фуллеренов с подложки. Все использован-

фуллерены, погруженные в алюминий, уменьшают

ные в модели эмпирические параметры были взяты

свой объем. Изначально при 1000 K радиус фулле-

из расчетов ab initio. Модель обеспечивает хорошее

рена R_f составлял 5.223Å (R = 3.705Å), а после

согласие с полученными в настоящей работе резуль-

погружения в расплав R_f = 5.207Å (R = 3.689Å).

татами расчетов методом МД и известными из ли-

Расчет по формуле (8) дает значение объемного мо-

тературы экспериментальными данными.

дуля K = 162 ГПа, которое отличается от результа-

Потенциал (1) был применен для исследования

тов DFT-вычислений [43,57] (868 ГПа и 874 ГПа со-

капиллярных явлений на границах между распла-

ответственно). Кроме того, наблюдается существен-

вом алюминия и фуллеренами. Рассчитанное при

ное расхождение результатов с данными работы [58],

температуре 1000 К значение давления Лапласа,

где с использованием различных эмпирических по-

действующего на фуллерены, составило 1.49 ГПа,

тенциалов для углерода были получены значения в

а коэффициент поверхностного натяжения γ

диапазоне от 370 ГПа до 694 ГПа. Указанные рас-

= 0.39 Дж/м². Это значение больше полученного

хождения лишь отчасти можно объяснить выбран-

нами ранее значения коэффициента поверхностного

ным в настоящей работе способом вычисления объ-

натяжения межфазной границы между алюминием

ема фуллерена. Использование в расчетах по фор-

и графитом [55].

мулам (5) и (8) в качестве радиуса фуллерена сред-

Поскольку расплав алюминия слабо смачивает

него расстояния от центра масс до ядер атомов дает

поверхность фуллерена, молекулы, распределенные

значение объемного модуля K_T = 324 ГПа. Это зна-

в расплаве, склонны к образованию связанных клас-

чение намного меньше величины 674 ГПа, получен-

теров. В настоящей работе коагуляция наблюдалась

ной с использованием потенциала Терсоффа в рабо-

в ходе прямого МД-моделирования, проводимого на

те [58]. Расхождение объясняется тем фактом, что в

временном интервале 10 нс. Характерное время про-

представленном исследовании использовалась ори-

цесса коагуляции на указанном временном интерва-

гинальная параметризация потенциала Терсоффа,

ле составляет τ1/2 ≈ 18.0 нс. Это намного больше

предложенная в работах [35, 36], в то время как ав-

времени релаксации диффузионного распределения

торы [58] модифицировали параметры потенциала

фуллеренов по объему, t_D ≈ 0.29 нс. Следователь-

с целью лучшего описания доступных эксперимен-

но, можно разделить «быстрый» процесс диффузи-

тальных данных. Несколько лучше рассчитанная в

онного распределения и «медленный» процесс коа-

настоящей работе величина K_T = 324 ГПа согла-

гуляции для приближенного анализа метастабиль-

суется со значением 370 ГПа из работы [58], полу-

ной системы с использованием методов равновесной

ченным с использованием потенциала [59]. Заметим,

статистической физики.

что в нашей работе объемный модуль был вычислен

при 1000 K, в то время как в [58] расчеты выполня-

RDF центров масс фуллеренов имеет ярко вы-

лись при нулевой температуре. Предполагая моно-

раженные максимумы и минимумы, что указыва-

тонное убывание величин упругих модулей с ростом

ет на сильное взаимодействие между молекулами.

температуры, мы считаем указанное соответствие

Взаимодействие носит капиллярный характер и осу-

ществляется в результате перекрывания полей упру-

разумным. Причины же расхождений полученных

результатов с данными DFT-моделирования могут

гих напряжений в расплаве алюминия.

быть связаны как с приближенным описанием меж-

Анализ кинетики коагуляции подтверждает на-

атомных взаимодействий при использовании потен-

личие сильного капиллярного взаимодействия меж-

циала [35], так и с приближенным учетом корреля-

ду фуллеренами. Применение модели Смолуховско-

ции электронов в DFT.

го ведет к завышению константы скорости коагу-

ляции примерно на два порядка величины по срав-

нению с рассчитанной при моделировании методом

4. ВЫВОДЫ

МД. Следовательно, процесс коагуляции имеет ак-

тивационный барьер, возникновение которого мож-

Исследовано взаимодействие поверхностей алю-

но объяснить наличием дальнодействующего потен-

миния в жидком и твердом состояниях с фуллерена-

циала межмолекулярного взаимодействия.

101

В. В. Решетняк, О. Б. Решетняк, А. В. Аборкин, А. В. Филиппов

ЖЭТФ, том 161, вып. 1, 2022

Рассчитаны значения КТР и объемного моду-

F. A. Khalid, O. Beffort, U. E. Klotz et al., Acta

ля фуллеренов. В интервале температур от 100 до

Mater. 51, 4575 (2003).

2000 K полученное значение КТР постоянно и рав-

A. I. Korobov, A. I. Kokshaiskii, V. M. Prokhorov et

но 1.15 · 10-5 K-1. Этот результат расходится с дан-

al., Phys. Solid State 58, 2472 (2016).

ными [56], где отмечена сильная зависимость β(T )

J. Shin, K. Choi, S. Shiko et al., Composites Part B:

при T ≤ 300 K, а при T ≥ 300 K — медленное из-

Eng. 77, 194 (2015).

менение значения β от 1.8 · 10-5 до 2.0 · 10-5 K-1.

Расхождения, полученные при 100 ≤ T ≤ 300 K, свя-

A. J. Maxwell, P. A. Brühwiler, S. Andersson et al.,

заны с использованием в настоящей работе потен-

Phys. Rev. B 52, R5546 (1995).

циала Терсоффа, который не воспроизводит зави-

A. J. Maxwell, P. A. Brühwiler, D. Arvanitis et al.,

симость β(T) при низких температурах. Причиной

Phys. Rev. B 57, 7312 (1998).

расхождения при T ≥ 300 K является выбор мето-

да расчета объема фуллерена. В настоящей работе

10.

M. K.-J. Johansson, A. J. Maxwell, S. M. Gray et al.,

объем вычислялся с использованием радиуса кон-

Phys. Rev. B 54, 13472 (1996).

тактной поверхности R_f , учитывая конечные значе-

11.

M. K.-J. Johansson, A. J. Maxwell, S. M. Gray et al.,

ния атомных размеров, в то время как в [56] было

Surface Sci. 397, 314 (1998).

использовано среднее расстояние R от центра масс

12.

R. Fasel, P. Aebi, R. G. Agostino et al., Phys. Rev.

фуллерена до атомов углерода. Пересчет КТР с ис-

Lett. 76, 4733 (1996).

пользованием R в качестве радиуса фуллерена дает

значение β = 1.63 · 10-5 K-1, которое согласуется с

13.

A. V. Hamza, J. Dykes, W. D. Mosley et al., Surface

данными из [56].

Sci. 318, 368 (1994).

Вычисленное значение объемного модуля фулле-

14.

M. Stengel, A. De Vita, and A. Baldereschi, Phys.

рена K_T = 162 ГПа. Это значение заметно отличает-

Rev. Lett. 91, 166101 (2003).

ся от имеющихся в литературе данных [43,57,58,60],

что лишь отчасти объясняется зависимостью объ-

15.

X. Q. Shi, M. A. Van Hove, and R. Q. Zhang, Phys.

Rev. B 85, 075421 (2012).

емного модуля от выбора способа расчета объема

наночастицы. Пересчет с использованием в качестве

16.

D. A. King, Surface Sci. 47, 384 (1975).

радиуса фуллерена среднего расстояния от атомов

17.

E. I. Altman and R. J. Colton, Surface Sci. 295, 13

углерода до центра масс молекулы дает величину

(1993).

324 ГПа. Анализ литературных данных [58] свиде-

тельствует также о чувствительности результатов

18.

H.-I. Li, K. Pussi, K. J. Hanna et al., Phys. Rev. Lett.

вычислений упругих модулей к выбору модели вза-

103, 056101 (2009).

имодействия для атомов углерода.

19.

S. Modesti, J. K. Gimzewski, and R. R. Schlittler,

Surface Sci. 331, 1129 (1995).

Финансирование. Работа выполнена при фи-

нансовой поддержке Министерства науки и высше-

20.

C. J. Villagomez, I. L. Garzon, and L. O. Paz-Borbón,

Comput. Mater. Sci. 171, 109208 (2020).

го образования Российской Федерации (соглашение

с ОИВТ РАН №075-15-2020-785).

21.

J. Hutter, M. Iannuzzi, F. Schiffmann et al., Wiley

Interdisciplinary Reviews: Comput. Molecular Sci. 4,

15 (2014).

ЛИТЕРАТУРА

22.

J. VandeVondele, M. Krack, F. Mohamed et al.,

Comput. Phys. Commun. 167, 103 (2005).

1. V. Chak, H. Chattopadhyay, and T. L. Dora, J.

Manuf. Proces. 56, 1059 (2020).

23.

G. Lippert, J. Hutter, and M. Parrinello, Theor.

Chem. Accounts 103, 124 (1999).

2. Y. Hu, O. A. Shenderova, Z. Hu et al., Soft Matter

24.

M. F. Peintinger, D. V. Oliveira, and T. Bredow, J.

3, 1099 (2007).

Comput. Chem. 34, 451 (2013).

3. S. R. Bakshi, D. Lahiri, and A. Agarwal, Int. Mater.

25.

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.

Rev. 55, 41 (2010).

Lett. 77, 3865 (1996).

4. H. G. P. Kumar and M. A. Xavior, Proc. Eng. 97,

26.

R. W. G. Wyckoff, Crystal Structures, Vol.

1033 (2014).

Interscience Publ. (1963).

102

ЖЭТФ, том 161, вып. 1, 2022

Межатомное взаимодействие на границе...

27.

H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13,

45.

Yu. A. Kvashnina, A. G. Kvashnin, L. A. Chernoza-

5188 (1976).

tonskii et al., Carbon 115, 546 (2017).

28.

J. Schöchlin, K. P. Bohnen, and K. M. Ho, Surface

46.

C. A. Perottoni and J. A. H. da Jornada, Phys. Rev.

Sci. 324, 113 (1995).

B 65, 224208 (2002).

29.

P. W. M. Jacobs, Yu. F. Zhukovskii, Yu. Mastrikov

47.

M. S. Dresselhaus, G. Dresselhaus, and P. C. Eklund,

et al., Comput. Model. New Technol. 6, 7 (2002).

Science of Fullerenes and Carbon Nanotubes: their

Properties and Applications, Elsevier (1996).

30.

C. Fiolhais, L. M. Almeida, and C. Henriques, Progr.

Surf. Sci. 74, 209 (2003).

48.

T. Liang, T.-R. Shan, Y.-T. Cheng et al., Mater. Sci.

Eng. R: Reports 74, 255 (2013).

31.

W. R. Tyson and W. A. Miller, Surface Sci. 62, 267

(1977).

49.

T. Liang, Y. K. Shin, Y.-T. Cheng et al., Ann. Rev.

Mater. Res. 43, 109 (2013).

32.

A. Kiejna, Phys. Rev. B 68, 235405 (2003).

33.

http://lammps.sandia.gov.

50.

I. A. Evdokimov, R. R. Khayrullin, R. Kh. Bagramov

et al., Russ. J. Non-Ferrous Metals 62, 132 (2021).

34.

S. Plimpton, J. Comput. Phys. 117, 1 (1995).

51.

A. V. Aborkin, A. I. Elkin, V. V. Reshetniak et al.,

35.

J. Tersoff, Phys. Rev. B 39, 5566 (1989).

J. Alloys Compounds 872, 159593 (2021).

36.

J. Tersoff, Phys. Rev. Lett. 61, 2879 (1988).

52.

L. Botto, E. P. Lewandowski, M. Cavallaro et al., Soft

Matter 8, 9957 (2012).

37.

H. W. Sheng, M. J. Kramer, A. Cadien et al., Phys.

Rev. B 83, 134118 (2011).

53.

M. V. Smoluchowski, Z. Phys. Chem. 92, 129 (1918).

38.

H. Heinz, R. A. Vaia, B. L. Farmer et al., J. Phys.

54.

B. V. Derjaguin, N. V. Churaev, and V. M. Muller,

Chem. C 112, 17281 (2008).

Surface Forces, Springer, Boston, MA (1987).

39.

L. S. Smith and L. L. Lee, J. Chem. Phys. 71, 4085

55.

V. V. Reshetniak and A. V. Aborkin, JETP 130, 214

(1979).

(2020).

40.

L. D. Landau and E. M. Lifshitz, Statistical Physics,

56.

Y.-K. Kwon, S. Berber, and D. Tománek, Phys. Rev.

Vol. 5, Elsevier Sci. (2013).

Lett. 92, 015901 (2004).

41.

W. Shinoda, M. Shiga, and M. Mikami, Phys. Rev.

57.

R. Peón-Escalante, C. Villanueva, R. Quintal et al.,

B 69, 134103 (2004).

Comput. Mater. Sci. 83, 120 (2014).

42.

A. P. Thompson, S. J. Plimpton, and W. Mattson, J.

Chem. Phys. 131, 154107 (2009).

58.

N. Kaur, S. Gupta, V. K. Jindal, and K. Dharamvir,

Carbon 48, 744 (2010).

43.

A. Khabibrakhmanov and P. Sorokin, Carbon 160,

228 (2020).

59.

D. W. Brenner, Phys. Rev. B 42, 9458 (1990).

44.

N. P. Kobelev, R. K. Nikolaev, Ya. M. Soifer et al.,

60.

R. S. Ruoff and A. L. Ruoff, Appl. Phys. Lett. 59,

Chem. Phys. Lett. 276, 263 (1997).

1553 (1991).

103