Моделирование взаимодействия нейтральных металлических нанокластеров при соударении с металлической поверхностью Текст научной статьи по специальности «Физика»

Математическое моделирование

УДК 531.395, 539.198, 539.231, 519.6

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

Б. Батгэрэл, Э. Г. Никонов, И. В. Пузынин

Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д.6, г.Дубна, Московская область, 141980, Россия

В работе приводятся результаты исследований методами классической молекулярной динамики процессов взаимодействия нейтральных металлических нанокластеров при соударении с металлической поверхностью. Исследована зависимость характерных размеров структуры образованного в результате соударения поверхностного слоя от размера, энергии соударения и частоты импульсного источника нанокластеров. В результате, численно определена функциональная зависимость глубины проникновения атомов кластера в материал мишени и толщины осаждаемого слоя от количества атомов в налетающих кластерах и частоты импульсного источника. Обнаружено также, что толщина осаждаемого слоя, в отличие от глубины проникновения, перестаёт зависеть от числа атомов в налетающих кластерах N, частоты импульсного источника ш и энергии налетающих кластеров Е при увеличении N, ш и Е. При этом осаждаемый слой становится неоднородным по толщине и принимает характерную воронкообразную форму. Показано, что существует зависимость характеристик различных энергетических режимов (soft landing, droplet spreading и implantation) от числа атомов в налетающих кластерах. Исследованные проблемы могут представлять интерес для развития технологий получения наноматериалов с новыми физическими и химическими свойствами.

Ключевые слова: молекулярная динамика, нанокластеры, моделирование, метод Верле, ударное взаимодействие.

1. Введение

В связи с развитием нанотехнологий увеличился интерес к моделированию процессов, происходящих при облучении материалов пучками высокоэнергетических наночастиц, в результате которых могут образовываться наноматериалы с новыми физическими, химическими и структурными свойствами [1—5]. Одним из распространённых в настоящее время подходов к численному моделированию этих процессов являются различные методы молекулярной динамики [5].

В работах [6-8] представлены результаты моделирования для металлов процессов модификации поверхности при облучении кластерами частиц с энергиями столкновения от (< 0.1 эВ/атом) до (> 3 эВ/атом). В работе [9] приведены результаты численных исследований, в рамках молекулярно-динамического моделирования, столкновения металлических нанокластеров между собой и с металлической поверхностью в широком интервале энергий столкновения. Обнаружен новый механизм образования связанного состояния в системах «кластер - кластер» и «кластер - поверхность» без плавления на контакте. Проведён также физический анализ этого процесса.

Целью данной работы является исследование методами классической молекулярной динамики процессов взаимодействия нейтральных металлических нано-кластеров при соударении с металлической поверхностью. В работе также исследована зависимость характерных размеров структуры образованного в результате

Статья поступила в редакцию 4 июня 2013 г. Работа выполнена при поддержке гранта РФФИ 12-01-00396а.

соударения поверхностного слоя от размера, энергии соударения и частоты импульсного источника нанокластеров. В результате численно определена функциональная зависимость глубины проникновения атомов кластера в материал мишени и толщины осаждаемого слоя от количества атомов в налетающих кластерах, частоты импульсного источника. Обнаружено также, что толщина осаждаемого слоя, в отличие от глубины проникновения, перестаёт зависеть от числа атомов в налетающих кластерах N, частоты импульсного источника ш и энергии налетающих кластеров Е при увеличении N, ш и Е. При этом осаждаемый слой становится неоднородным по толщине и принимает характерную воронкообразную форму. Показано, что существует зависимость характеристик различных энергетических режимов^оИ landing, droplet spreading и implantation) от числа атомов в налетающих кластерах.

2. Методы молекулярной динамики

Молекулярно-динамическое моделирование основано на представлении объектов моделирования в виде системы взаимодействующих частиц (атомов или молекул). Эволюция системы происходит в результате движения упомянутых выше частиц. Координаты частиц в каждый последующий момент времени вычисляются посредством интегрирования уравнений движения, в которые входят потенциалы взаимодействия частиц между собой и внешней средой.

В зависимости от типа уравнений движения молекулярная динамика подразделяется на следующие независимые системы моделирования: классическая молекулярная динамика, квантовая молекулярная динамика, релятивистская молекулярная динамика, релятивистская квантовая молекулярная динамика. Характерные значения энергий, расстояний и временных промежутков, в процессах столкновения нейтральных металлических нанокластеров с металлической поверхностью позволяют проводить моделирование описанной системы в рамках классической молекулярной динамики.

В основе методов классической молекулярной динамики лежит модельное представление о многочастичной системе, в которой все частицы (атомы или молекулы) представляют собой материальные точки. Поведение отдельной частицы описывается классическими уравнениями движения Ньютона [10], которые могут быть записаны в следующем виде.

^ = ^ (1)

Здесь % — номер частицы (1 ^ % ^ N), N — полное число частиц, mi — масса

частицы, fi — равнодействующая всех сил, действующих на частицу, имеющая следующее представление.

? _ ди(П,--,гЪ) , ?, .

Ь = ^ + !г , (2)

->ех

где и — потенциал взаимодействия между частицами, Д — сила, обусловленная внешними полями.

Потенциал взаимодействия между частицами может быть задан в виде силового поля [11].

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

и = ^ Ubond + ^ Uangle + ^ Udihed + ^ Uinv + У^ U pair + — body — body +

^ ^ U^ersof f + ^ ^ Umetal.

Здесь иЬопа — потенциал, описывающий непосредственное прямое взаимодействие между атомами (частицами); иапдIе — потенциал, описывающий непрямое взаимодействие между атомами (частицами); нел — потенциал, описывающий взаимодействие между атомами (частицами), возникающее в результате возникновения торсионных сил внутри молекулы; Щпг, — потенциал, описывающий взаимодействие между атомами (частицами), возникающее из-за особенностей геометрического расположения трёх атомов вокруг некоторого центрального внутри молекулы; ира1 г, и3-ъо а у — многочастичные потенциалы; и4-ъо ау, итегзо// — потенциалы, используемые для описания многокомпонентных систем; итега1 — металлический потенциал.

В данной работе для моделирования взаимодействия пучком нанокластеров с металлической поверхностью использовался потенциал Саттона -Чена.

1 N N N N

итеЛаI = 2 X) X УИ(ГИ) + X Р(Р*)' Р* = X Р^(гЧ)'

1=1 j=i 1=1 Ъ=1, i=j

(3)

где = € ^^^ — потенциал парного взаимодействия, рц(тц) = —

вклад в плотность заряда г-го атома от ^-го атома, Е(р^) = —с€^/~р1 — функция «погружения», которая представляет собой энергию, необходимую для перемещения г-го атома в электронное облако.

Значения величин параметров а,с,€,т и п для потенциала Саттона -Чена приведены в табл. 1.

Таблица 1

Значения параметров для потенциала Саттона—Чена

Ме1а1 а(А) с €(еУ) т п

А1 4.05 16.399 0.0331 6 7

Си 3.61 39.755 0.0124 6 9

В реальных экспериментах температура системы поддерживается за счёт энергообмена с внешней средой. Детальный учёт взаимодействия частицы с внешней средой часто невозможен. Для учёта эффектов энергообмена с внешней средой используются специальные алгоритмы - термостаты. В молекулярной динамике температура молекулярной системы вводится через удельное среднее значение кинетической энергии. Выражение для средней кинетической энергии системы имеет вид:

Е =

2 N

Е т^2

¿=1

N

(4)

где mi — масса частицы, Ьг — скорость частицы, N — полное число частиц.

Из статистической физики известно, что кинетическая энергия системы и её температура связаны следующим соотношением:

Е =

3кТ

Т~'

(5)

к — постоянная Больцмана.

Из (4) и (5) получим мгновенное значение температуры:

Т =

N

Е

г=1_

3Мк '

(6)

Далее, проведя усреднение по времени, получим значение температуры системы:

t0+T

N

to i=1

Использование термостата особенно важно на этапе релаксации системы. В случае установившегося термодинамического равновесия температура термостата и средняя температура молекулярной системы должны совпадать. Энергии подсистем обычно много меньше энергии термостата, это является условием приближенного термодинамического равновесия. При моделировании методами молекулярной динамики обычно фиксируют температуру термостата. Температура молекулярной системы может при этом меняться вследствие различных причин. Например, из-за конечного шага интегрирования частица может оказаться в классически запрещённой области. Это приведёт к резкому скачку энергии, а затем и температуры. Наиболее часто встречающимися моделями термостатов являются коллизионный термостат, основанный на столкновительной динамике, термостат Берендсена, использующий в уравнениях движения знакопеременное нелинейное трение, термостат Эванса, Ланжевена, Андерсена, Nose - Hoover термостат и другие. В частности, термостат Эванса имеет следующее определение.

d^l m dvj(t) Ш m m

~dT =Vi(t), ~dT = -w - x(t)Vi(t), (R)

di l^mivi(t)

В данной работе для учёта эффектов теплообмена системы «кластер-мишень» с внешней средой использовался термостат Берендсена, который реализуется при помощи алгоритма (8). При этом для обеспечения слабой связи с внешней средой выражение для коэффициента пересчёта скоростей \(t) имеет следующий вид.

x(t) =

At i ТР т+

1 +--_ 1

+ ТВ \Т(t)

(9)

характерные значения постоянной времени находятся в интервале тв G [0.5,2] псек.

Для обеспечения термодинамического равновесия при моделировании сложных многокомпонентных систем, как правило, также фиксируется давление в системе. Для решения данной задачи используются так называемые баростаты. Наиболее часто используемыми моделями баростатов являются баростат Беренд-сена, баростат Ланжевена и Nose - Hoover баростат.

Для интегрирования системы уравнений (1) необходимо задать начальные координаты и импульсы частиц и граничные условия. При этом для корректного интегрирования уравнений движения молекулярной динамики необходимо учитывать следующее.

Длина траектории в молекулярной динамике равняется шагу интегрирования, умноженному на число произведённых шагов. Выбор длины траектории в значительной степени связан с понятием эргодичности траектории. В молекулярной динамике обычно имеют дело со средними величинами вдоль траектории (или со средними по времени). В эксперименте обычно имеют дело с величинами средними по ансамблю. Для того чтобы сравнение статистических характеристик системы с результатами молекулярно-динамических расчётов было корректным, необходимо, чтобы траектория обладала достаточно хорошими эргодическими свойствами [10]. Реально это означает, что за время интегрирования система должна много раз побывать во всех значимых областях конфигурационного пространства, т.е. время интегрирования должно быть значительно больше т - времени

2

преодоления барьеров.

г - 3м ехр ^' (10)

где М — количество торсионных углов в молекуле, и — значение энергетического барьера, к — постоянная Больцмана, Т — температура.

Используя (10), можно оценить минимальную длину траектории, которая должна быть значительно больше, чем время, необходимое для преодоления каждого из энергетических барьеров.

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

fi = -Уг X U(j'ij). (11)

Затем рассчитываются новые координаты частиц, из которых определяются равнодействующие силы:

ri(t + At) = fi(t) + Vi(t)At + ai^At2. (12)

Здесь a — ускорение, a(t + At) = ^.

Далее определяются скорости атомов:

ait + At) + ait) . v(t + At) = v(t) +-i-2-^At. (13)

Одной из наиболее существенных проблем процедуры интегрирования является выбор шага. При большом шаге погрешности интегрирования могут быть значительными, что приведёт к нестабильности траектории. При малом шаге существенно увеличивается время расчёта. В связи с этим шаг интегрирования был выбран равным At = 10- 15сек = 1фс. В уравнениях движения, описывающих изменения по различным степеням свободы, временные характеристики существенно отличаются друг от друга. Для достаточно точного вычисления решения по быстрым и медленным переменным шаги интегрирования по ним могут различаться. По быстрым переменным может быть выбран значительно больший шаг [10]. В методе Верле шаг интегрирования берётся единым. Начальные скорости атомов выбираются с помощью генератора случайных чисел в соответствии с распределением Максвелла при заданной температуре.

Наиболее используемыми модификациями метода Верле являются следующие.

Метод Velocity Verlet.

В рамках метода Velocity Verlet интегрирование уравнений движения производится по следующей схеме:

— в начале каждого шага задаются или рассчитываются на предыдущем шаге величины t : r(t), v(t), f (t);

— затем вычисляются значения скоростей частиц в момент времени t + ^ и координат нового местонахождения частиц

/ At \ . At f (t) . , А f At \

vit + —j = v(t) + A • ^, r(l + At) = r(t) + At • vit + —j ;

— далее пересчитываются силы, действующие на частицу в момент времени t + At,

f (t + At) = f (t);

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

( At\ At f (t + At) vif + At) = ^t +T) .

Метод Leapfrog Verlet.

Метод Leapfrog Verlet состоит из одного этапа и реализуется по следующей схеме:

— в начале каждого шага задаются или рассчитываются на предыдущем шаге величины t : r(t), f (t — At) и скорости v(t — 4 ) в момент времени t — 4 ;

— далее пересчитываются силы, действующие на частицу в момент времени t — At,

f (t) = f (t — At);

— далее вычисляются значения скоростей в момент времени t + 4

f (t)

(t + = v(t - + At • ^, r(t + At) = r(t) + At • v(t + A2pj

2 ) \ 2 ! m

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

(14)

v(t)=2

v(* - t) + v{* +A)

Для анализа результатов моделирования рассчитываются различные характеристики системы. Кинетическая энергия:

Температура:

где степень свободы:

l

N

Ekin(t) = míV¡(t).

=1

T (t) = — Ekin, ka

(15)

(16)

d 3N 3Nfrozen ^constraints 3 P,

(17)

для периодического граничного условия: р = 3, для непериодического граничного условия: р = 0.

3. Программное обеспечение для молекулярно-динамического программирования

В данной работе компьютерное моделирование проводилось с использованием пакета программ Las Palmeras Molecular Dynamics (LPMD) [12]. Данный пакет, в отличие от множества программных продуктов, предназначенных для молеку-лярно-динамического моделирования (DL_POLY [l1], Gromacs, Moldy, Moscito и др.), является наиболее подходящим для исследования проблем, сформулированных в данной работе.

Программа LPMD предназначена для молекулярно-динамического моделирования как равновесных, так неравновесных процессов, таких как распространение ударных волн, бомбардировка снарядами мишени, кластерных столкновений, распространение тепла и т.п.

Пакет реализован на языке программирования CH—+, имеет модульную архитектуру. Функциональность программной системы может быть расширена за счёт

подключения p1ugin-модулей, которые вызываются по мере необходимости во время запуска основной программы. Данный программный продукт содержит также набор модулей для визуализации процессов моделирования в реальном времени.

4. Результаты моделирования

Для моделирования использовалась следующая конфигурация системы «снаряд-мишень». На неподвижную мишень налетает пучок нанокластеров, сформированных из частиц, обладающих гранецентрированной структурой (рис. 1).

Мишень представляет собой металлический прямоугольный параллелепипед с

—

Рис. 1. Схема эксперимента

геометрическими размерами 54.15 х 108.3 х 108.3 А и состоит из 54000 атомов меди. Моделирование проводилось для наночастиц, состоящих из 13, 55 и 147 атомов меди с поперечными размерами 5.86 А, 12.17 А и 15.82 А соответственно. Рассматривались два варианта: облучение мишени одиночными наночастицами и облучение мишени последовательностью наночастиц при энергиях налетающих кластеров равных 0.1, 1.0, 10 эВ/атом. Непосредственные вычисления проводились в два этапа. Сначала перед непосредственным столкновением проводилась специальная процедура подготовки кластеров. Затем проводилось моделирование процессов столкновения сформированных нанокластеров с твёрдой поверхностью.

4.1. Процедура подготовки нанокластеров

Данная процедура необходима для учёта реальных условий образования пучков кластеров в экспериментах, поскольку структура реальных наноразмерных кластеров существенно отличается от гранецентрированной кристаллической решётки. Сначала из достаточно большой гранецентрированной кристаллической решётки выбирается необходимое для получения оптимальной (как правило, ико-саэдрической) структуры нанокластера число (13, 55, 147,...) атомов, находящихся на наименьших расстояниях от центрального атома. В начальном состоянии наночастицы представляют собой частицы, имеющие структуру гранецентрированной кубической кристаллической решётки размером 14.44 A . Следующим шагом является термализация наночастиц без учёта граничных условий. Затем постепенно нагреваем получившийся кластер до температуры Т = 1800°К в течение 105 шагов по времени (шаг по времени равен 1 фс, температура плавления меди — Т = 1356°К), потом охлаждаем до температуры Т = 300°К и проводим процедуру термализации. Данная процедура проводится также в течение 105 шагов с шагом 1 фс. Достижение кластером оптимальной структуры подтверждается установлением парных структур типа 1-5-5-5, характерных для икосаэдрических решёточных структур, что выявляется в рамках структурного анализа методом CAN (common neighbor analysis) [12]. Суть данного метода состоит в анализе взаимного расположения соседних атомов с использованием кристаллографических

индексов. Каждая пара атомов характеризуется наборами чисел четырёх индексов (г,'}, к, I). Здесь индекс % равен 1 для близких соседей и 2 для остальных, второй индекс '} равен числу общих близких соседей, третий индекс к равен числу бондов (ковалентных связей между атомами) среди ] общих близких соседей и четвёртый индекс I равен длине наибольших связанных бондов среди общих близких соседей. Распределение структур с различными кристаллографическими индексами, демонстрирующими преобладание икосаэдрических решёточных структур в результате процедуры подготовки нанокластеров, представлено в табл. 2 и 3.

— относительное количество кристаллических структур с определёнными индексами в процентах.

Таблица 2

Распределение кристаллических структур в первоначальных наночастицах с преобладанием гранецентрированных структур

1 .1 к 1 ^ам

1 1 0 0 0.15

1 2 0 0 7.22

1 2 1 1 7.22

1 3 1 1 23.42

1 4 2 1 62.00

Таблица 3

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

нанокластеров

1 1 к 1 ^•ам

1 0 0 0 1.03

1 1 0 0 1.47

1 2 0 0 6.45

1 2 1 1 6.30

1 3 0 0 2.79

1 3 1 1 18.48

1 3 2 2 5.13

1 4 1 1 0.29

1 4 2 1 28.01

1 4 2 2 27.27

1 4 3 3 1.47

1 5 4 4 1.17

1 5 5 5 0.15

4.2. Моделирование процессов столкновения

После образования нанокластеров выполнялось численное моделирование процессов столкновения с металлической поверхностью. Сначала моделировалось столкновение одиночных нанокластеров, затем столкновение пучка их трёх кластеров с поверхностью. Моделирование проводилось для трёх типов кластеров, состоящих из N = 13, 55 и 147 атомов для значений энергии Е налетающих нано-частиц 0.1, 1.0, и 10 эВ/атом. На рис. 2 и 3 приведены для примера результаты моделирования процессов соударения одиночного кластера для числа атомов в кластере N = 147 и при энергии Е =10 эВ/атом. В данном варианте расчётов структурные изменения мишени и кластеров наиболее очевидны по сравнению с вариантами расчётов при меньших значениях N и Е.

В результате моделирования процессов бомбардирования поверхности пучками наночастиц, состоящими из трёх кластеров, мы получаем следующую картину.

Рис. 3. Кластеры N = 147, Е =10 эВ/атом, t =1.2,1.3,1.5, 2, 3,4, 5,10 пс. Вид

спереди

Рис. 4 и 5 демонстрируют результаты моделирования для кластеров, содержащих N = 147 атомов для значений энергии Е налетающих наночастиц 10 эВ/атом.

После облучения поверхности мишени нанокластерами исследовалась зависимость глубины проникновения h составляющих частиц нанокластеров в материал мишени и толщины d поверхностного слоя, образованного в результате облучения, от величины нанокластера и энергии налетающих нанокластеров после термали-зации всей системы (рис. 6).

Исследования глубины проникновения нанокластеров в облучаемый материал проводилась для кластеров, содержащих N = 13, 55 и 147 атомов для значений энергии Е налетающих наночастиц 0.1, 1.0, и 10 эВ/атом, соответственно.

Рис. 4. Пучки кластеров N = 147, Е = 10 эВ/атом,

t = 1.3, 2, 4, 6.4,10,11.4,15, 20 пс. Вид сбоку

Рис. 5. Пучки кластеров N = 147, Е = 10 эВ/атом,

t = 1.3, 2,4, 6.4,10,11.4,15, 20 пс. Вид спереди

Рис. 6. d — толщина поверхностного слоя, h — глубина проникновения частиц нанокластеров в материал мишени

На рис. 7 приведены графики зависимости глубины проникновения нанокластеров в облучаемый материал h от времени для числа частиц в кластере N = 147 и энергии налетающих частиц Е = 10 эВ/атом. Данные графики наиболее наглядно демонстрируют основные особенности процессов проникновения атомов нанокластеров в материал мишени. На рис. 8 приведены графики зависимости глубины проникновения h и толщины поверхностного слоя d от энергии Е налетающих нанокластеров 0.1, 1.0, и 10 эВ/атом для кластеров, содержащих N = 13, 55

и 147 атомов и для случаев одиночных нанокластеров (Single, Т = то) и пучка из трёх нанокластеров (Т = 4, 5, 6 ps).

Рис. 7. Зависимость к глубины проникновения нанокластеров в облучаемый материал от времени для N = 147 и Е =10 эВ/атом

5. Заключение

Задача исследования зависимости плотности и температуры от энергии в системе кластер-поверхность далеко не тривиальная задача [6-8]. Тем не менее мо-лекулярно динамическое моделирование позволяет описывать такого типа задачи с высоким разрешением по пространственным и временным координатам. Для удобства проведения анализа результатов моделирования мишень разбивают на конечное количество слоев фиксированной толщиной ао вдоль движения налетающих кластеров (обычно координатная ось z). Все полученные данные (для координат положения частиц, скоростей и т.д.) анализируют в каждом из таких слоев (например, вычисление средней величины скорости и т.п.). Для каждого такого слоя рассчитывается изменение величин плотности и температуры на каждом временном шаге для мишени и снаряда. Анализ результатов моделирования позволяет сделать вывод о наличии трёх интервалов энергии столкновения кластеров с поверхностью, которые существенно отличаются по характеристикам основных процессов, происходящих при взаимодействии кластеров с материалом поверхности. Режимы соударения, соответствующие трём интервалам энергии, в литературе обычно называют следующим образом. При энергиях налетающих кластеров 0 < Е < 0.1эВ/атом происходит Soft landing, или налипание кластеров на поверхность. При энергиях Е порядка 1эВ/атом происходит Droplet spreading, или осаждение кластеров на поверхность с частичной деформацией налетающих кластеров. И наконец, при энергиях Е > 10эВ/атом происходит Implantation или глубокое проникновение кластеров в материал поверхности с образованием структуры поверхностного слоя. Процессы, происходящие при соударениях кластеров с поверхностью для значений энергий из описанных выше интервалов, характеризуются существенно разными распределениями плотности и температуры, а также величинами энергии осаждения.

12 Г

N=13 Ю - Single T = 4 ps T = 5 ps 8 - T = 6 ps

6 -

-s

4 2

4 5 6 E (eV/atom)

12 -,—

10 -

8 -

6 -

4 -

N=13 Single T = 4 ps T = 5 ps T = 6 ps

0

3 4 5 6 E (eV/atom)

0

0

2

3

7

8

9 10

0

2

7

8

9 10

Рис. 8. Зависимость глубины h проникновения нанокластеров и толщины d

поверхностного слоя от энергии

Soft landing.

Для энергии налетающих кластеров 0 < Е < 0.1эВ/атом поверхность мишени не разрушается и структура кластера остаётся очень близкой к первоначальному состоянию. Общим свойством процесса налипания является то, что кластер имеет локализованные распределения плотности и температуры и температура в данном процессе всегда ниже температуры плавления. В начале процесса соударения температура кластера во внешнем слое возрастает до своего максимального значения. Плотность кластера увеличивается примерно в два раза по сравнению с первоначальным значением аналогично результатам работ [6-8]. Как и ожидается, для твёрдого тела, не наблюдается видимых изменений в распределении плотности в материале мишени в случае облучения одиночными кластерами, независимо от размера кластеров. Аналогичная картина наблюдается и в случае облучения пучками кластеров. Глубина проникновения h при данных энергиях столкновения слабо зависит от числа атомов в налетающем кластере и практически равна нулю (рис. 7, 8). Напротив, толщина осаждаемого слоя d существенно отличается от нуля и растёт с ростом числа атомов в кластере и падает с ростом частоты импульсного источника пучка кластеров.

Droplet spreading.

При энергиях Е порядка 1эВ/атом зависимость температуры и плотности материала мишени становится неоднородной вдоль направления движения кластеров по сравнению с soft landing режимом и не восстанавливается до первоначальных значений температуры и плотности, которые были до соударения. При этом увеличивается глубина проникновения h и уменьшается толщина осаждаемого слоя d (см. рис. 7, 8). Кроме того, появляется явная зависимость глубины проникновения h от числа атомов в налетающем кластере.

Implantation.

При энергиях Е > 10эВ/атом становятся очевидными изменения структуры материала мишени после соударения с кластерами (см. рис. 2, 3). Причём очевидно также, что структурные изменения при облучении пучками кластеров являются значительно более существенными (см. рис. 4, 5). Аналогичные выводы можно сделать и в результате анализа зависимости температуры и плотности от времени. Глубина проникновения h при данных энергиях продолжает увеличиваться и явно зависит от числа атомов в налетающих кластерах и частоты импульсного источника в случае облучения мишени пучками кластеров. Толщина осаждаемого слоя d перестаёт увеличиваться с увеличением числа атомов в налетающих кластерах и частоты импульсного источника (см. рис. 7, 8). При этом осаждаемый слой становится неоднородным по толщине и принимает характерную воронкообразную форму (рис. 9).

Рис. 9. Центральное сечение материала мишени, перпендикулярное

фронтальной плоскости

Перспективы дальнейших исследований.

Как показали исследования, проведённые в данной работе, процессы, происходящие в системе снаряд-мишень, в случае облучения одиночными кластерами и в случае облучения пучками кластеров существенно отличаются по ряду параметров. Обнаружена зависимость глубины проникновения атомов кластера в материал мишени и толщины осаждаемого слоя от количества атомов в налетающих кластерах, частоты импульсного источника в случае облучения пучками кластеров. Обнаружено также, что толщина осаждаемого слоя, в отличие от глубины проникновения, перестаёт зависеть от числа атомов в налетающих кластерах N, частоты импульсного источника ш и энергии налетающих кластеров Е при увеличении N, ш и Е. Результаты моделирования показывают зависимость характеристик различных энергетических режимов (soft landing, droplet spreading и implantation) от числа атомов в налетающих кластерах. Таким образом, увеличивая число атомов в налетающих кластерах, можно переходить от одного режима к другому. Дальнейшие исследования позволят получить более полную и точную картину описанных выше эффектов. Кроме того, на наш взгляд, вызывает несомненный интерес исследование характеристик системы снаряд-мишень для различных параметров металлических потенциалов и выяснение условий переходов между различными режимами процессов столкновений. При этом увеличение

числа атомов в несколько раз в налетающих нанокластерах и на несколько порядков в материале мишени позволит производить реалистичные расчёты для наблюдаемых величин с учётом динамических и стохастических свойств молекулярно-динамических систем, связанных с локальной неустойчивостью траекторий и погрешностями численного интегрирования [13].

Литература

1. Макаров Г. Н. Экстремальные процессы в кластерах при столкновении с твердой поверхностью // Успехи физических наук. — 2006. — Т. 176, № 2. — С. 121174. [Makarov G. N. Extreme Processes in Clusters Impacting on a Solid Surface // UFN. — 2006. —Vol. 176, No 2. — P. 121-174. ]

2. Бойко В. И., Валяев А. Н., Погребняк А. Д. Модификация металлических материалов импульсными мощными пучками частиц // Успехи физических наук. — 1999. — Т. 169, № 11. — С. 1243-1271. [Bojko V.I., Valjaev A.N., Pogrebnjak A. D. Metal Modification by High-Power Pulsed Particle Beams // UFN. — 1999. — Vol. 169, No 11. — P. 1243-1271. ]

3. Диденко А. Н, Лигачев А. Е., Куракин И. Б. Воздействие пучков заряженных частиц на поверхность металлов и сплавов. — М.: Энергоатомиздат, 1987. — 183 с. [Didenko A.N., Ligachev A. E., Kurakin I. B. The Interaction of Charged Beams with the Metals and Alloys Surfaces. — Moscow: Energoatomizdat, 1978. — 183 p. ]

4. Действие излучения большой мощности на металлы / С. И. Анисимов, Я. А. Имас, Г. С. Романов, Ю. В. Ходыко. — М.: Наука, 1970. — С. 272. [Action of High-Power Radiation on Metals / S.I. Anisimov, Ya. A. Imas, G. S. Romanov, Yu. V. Khodyko. — Moscow: Nauka, 1970. ]

5. Быковский Ю. А., Неволин В. Н., Фоминский В. Ю. Ионная и лазерная имплантация металлических материалов. — М.: Атомиздат, 1991. — С. 235. [Bikovskiy Yu. A., Nevolin V. N., Fominskiy V. Yu. Ion and Laser Implantation of Metallic Materials. — Energoatomizdat, Moscow, 1991. ]

6. MD Simulation of Cluster-Surface Impacts for Metallic Phases: Soft Landing, Droplet Spreading and Implantation / K. Kholmurodov, I. Puzynin, W. Smith et al. // Computer Physics Communications. — 2001. — Vol. 141. — Pp. 1-16.

7. Molecular Dynamics Simulation of Cluster-Beam-Surface Impact Processes for Metallic Phases / K. Kholmurodov, I. Puzynin, W. Smith et al. // Journal of Computational Methods in Sciences and Engineering. — 2002. — Vol. 2, No 1s-2s. — Pp. 141-147.

8. Методы молекулярной динамики для моделирования физических и биологических процессов / Х. Т. Холмуродов, М. В. Алтайский, И. В. Пузы-нин и др. // Физика элементарных частиц и атомного ядра. — 2003. — Т. 34, № 2. — С. 472-515. [Molecular Dynamics Methods for Simulation of Physical and Biological Processes / H. T. Holmurodov, M. V. Altajskij, I. V. Puzynin et al.// Fizika elementarnykh chastic i atomnogo yadra. — 2003. — Vol. 34, No 2. — P. 472-515. ]

9. Головнев И. Ф., Головнева Е. И., Фомин В. М. Молекулярно-динамическое исследование столкновения нанокластеров друг с другом и с подложкой // Физическая мезомеханика. — 2007. — Т. 10, № 2. — С. 5-13. [Golovnev I.F., Golovneva E. I., Fomin V.M. Molecular Dynamics Investigation of Cluster Collisions with Each Other and with a Substrate // Fizicheskaja mezomehanika. — 2007. — Vol. 10, № 2. — P. 5-13. ]

10. Шайтан К. В., Терёшкина К. Б. Молекулярная динамика белков и пептидов. Методическое пособие. — 1999. — http://www.moldyn.ru/library/manual/ index.htm. [Shaitan K. V., Terjoshkina K. B. Molecular Dynamics of Proteins and Peptides. Methodical textbook. — 1999. ]

11. Smith W, Yong C. W., Rodger P. M. DL_POLY: Application to Molecular Simulation // Molecular Simulation. — 2002. — Vol. 28, No 5. — Pp. 385-471.

12. Las Palmeras Molecular Dynamics: Flexible and Modular Molecular Dynamics / S. Davis, C. Loyola, F. Gonzalez, J. Peralta // Computer Physics Communications. — 2010. — Vol. 181, No 12. — Pp. 2126-2139.

13. Норман Г. Э., Стегайлов В. В. Стохастическая теория метода классической молекулярной динамики // Математическое моделирование. — 2012. — Т. 24, № 6. — С. 3-44. [Norman G. Je., Stegajlov V. V. Stochastic Theory of the Classical Molecular Dynamics Method // Matem. modelirovanie. — 2012. — Vol. 24, No 6. — P. 3-44. ]

UDC 531.395, 539.198, 539.231, 519.6

Simulation of Impact Interaction of Uncharged Metallic Nanoclusters with Metallic Surface

B. Batgerel, E. G. Nikonov, I. V. Puzynin

Laboratory of Information Technologies Joint Institute for Nuclear Research Joliot-Curie, 6, 141980 Dubna, Moscow region, Russia

Investigation results of impact interaction of uncharged metallic nanocluster with metallic surface are presented. Simulation and investigations of impact processes are fulfilled by molecular dynamics methods and suitable software. Characteristic dimensions of surface layer produced by impact as a functions of cluster size, colliding energy and frequency of impulsive nanoclusters source have been analysed. It was found out that penetration depth and deposited layer thickness depend on number of particles in colliding nanoclusters and frequency of impulsive nanoclusters source. It was also discovered that deposited layer thickness in contrast to penetration depth ceases depending on number of particles in colliding nanoclusters N, and frequency of impulsive nanoclusters source w and colliding energy E with increasing of N, w and E. And at the same time deposited layer becomes heterogeneous in thickness and gets a funnel-shaped form. It is shown that realization of one of the choice of nanoclusters surface interaction (soft landing, droplet spreading and implantation) should be controlled by means of changing of both nanoclusters beam energy and number of atoms in clusters. Investigation results should be of interest in various fields of technologies developing nanomaterial with new physical and chemical properties.

Key words and phrases: molecular dynamics, simulation, Verlet method, impact interaction.