Редактирование: Моделирование поверхностной диффузии кремния методом молекулярной динамики

== Введение ==
В рассматриваемой работе описывается компьютерная модель,
позволяющая рассчитывать поведение кристалла кремния, имеющего
атомную решетку структуры алмаза. Для построения модели
использовался описанный ниже метод молекулярной динамики.
Приведены результаты численных экспериментов: нагрев кристалла с
измерением кинетической и полной энергий, нагрев кристалла с
свободной поверхностью 1,0,0, поведение атома на поверхности
1,1,1 и моделирование диффузии атомов на поверхности 1,1,1.

== Описание метода молекулярной динамики ==
В основе применения метода молекулярной динамики для численных
экспериментов лежит интегрирование уравнений движения частиц, –
атомов или молекул. При этом могут учитываться как потенциальные
силы, действующие между отдельными частицами, так и внешние
взаимодействия, например, для обеспечения нагрева и охлаждения
рассматриваемой среды.

== Задание свойств атомов кремния ==
Силы взаимодействия между атомами кремния в кристалле весьма
сложны, что обусловлено структурой электронных оболочек атомов.
Поскольку полная энергия изолированной системы должна сохраняться,
силы взаимодействия носят потенциальный характер. Т.е. сила,
действующая на один атом, может быть представлена как градиент от
некоторого потенциала – функции, зависящей от относительного
положения частиц. Наиболее простая разновидность потенциала
взаимодействия предполагает разбиение его на сумму слагаемых,
каждое из которых зависит только от расстояния между какими-либо
двумя атомами. Такую структуру имеют потенциалы [[Потенциал Леннарда-Джонса|Леннарда-Джонса]],
[[Потенциал Ми|Ми]], [[Потенциал Морзе|Морзе]].
Для них сила взаимодействия между двумя атомами может
быть вычислена с помощью дифференцирования потенциала по
скалярному аргументу – расстоянию между этими двумя атомами.
Подобные потенциалы не могут описать поведение атомов в
рассматриваемом кристалле. В описываемой компьютерной модели
используется [[Потенциалы Терсоффа, Бреннера|эмпирический потенциал взаимодействия Терсофа]] для
атомов кремния [1], [2]. Потенциал Терсоффа относится к многочастичным
потенциалам [3]. Так же, как вышеупомянутые потенциалы, многочастичные
потенциалы представляются в виде суммы слагаемых по всем парам
частиц, но каждое слагаемое зависит и от положения других частиц.

== Описание компьютерной модели ==
Компьютерный эксперимент производится посредством вычисления
радиус векторов и векторов скорости частиц в зависимости от
времени. Интегрирование ведется методом центральных разностей.
Метод состоит в том, что координаты и силы вычисляются во
временных точках, разделенных интервалами, равными шагу
интегрирования, а скорости вычисляются во временных точках,
находящихся в серединах вышеупомянутых интервалов:

<math>
    \bf{\underline{v}} (t + \tau / 2) = \underline{v} (t - \tau / 2) + \underline{w} (t) \tau
</math>

<math>
    \bf{\underline{r}} (t + \tau) = \underline{r} (t) + \underline{v} (t + \tau / 2) \tau
</math>

где <math>\tau</math> – шаг интегрирования. Ускорение <math>\bf{\underline{w}}(t)</math>
вычисляется через силу в тот же момент времени. В начальный
момент времени задается отсчетная конфигурация. Это может быть
кристалл, представляющий собой пластину с поверхностью 1,0,0 или
1,1,1. Конфигурация представляет собой идеальную решетку структуры
алмаза. Межатомное расстояние вычисляется аналитически, с помощью
минимизации функции потенциальной энергии ([[Потенциалы Терсоффа, Бреннера|потенциал Терсоффа]]). Для идеальной
решетки все слагаемые <math>V_{ij}</math> равны между собой. Кроме того,
для всех углов <math>\theta_{ijk} \quad \cos \theta_{ijk} = 1/3</math>. Это позволяет
вычислить коэффициент <math>b_{ij}</math>, как
постоянную величину, после чего выражение для <math>V_{ij}</math> минимизируется по
скалярному аргументу <math>r_{ij}</math>. Скорости в начальный момент
времени задаются как случайные величины в некотором диапазоне.
Величина диапазона зависит от начальной температуры, которую мы
хотим сообщить кристаллу. Функция распределения по скоростям в
начальный момент времени особого значения не имеет, поскольку
через промежуток времени порядка нескольких минимальных периодов
колебаний решетки устанавливается максвелловское распределение
скоростей. Перед интегрированием задается область пространства
– бокс, имеющий форму параллелепипеда, в котором будет проводиться
интегрирование. Число частиц и координаты бокса остаются
неизменными. Основной цикл программы – цикл интегрирования по
времени. Внутри каждого цикла вызывается функция, осуществляющая
один шаг интегрирования по времени. Шаг интегрирования выбирается
экспериментально, из соображений устойчивости вычислительной
схемы. Необходимым условием является сохранение полной энергии при
отсутствии внешних сил. В проводимых экспериментах использовалась
величина шага интегрирования, равная <math>\tau = 0.025</math> в единицах измерения,
наиболее удобных для использования в вычислительной схеме. Так,
координаты измеряются в ангстремах, единица времени составляет

<math>\mbox{\AA} \sqrt{\frac{m}{\mbox{eV}}} = 5.39523 \cdot 10^{-14}</math> с

где <math>m = 28.0855 \cdot 1.66054 \cdot 10^{-27}</math> кг –
масса атома кремния, <math>\mbox{\AA} = 10^{-10}</math> м – ангстрем,
<math>\mbox{eV} = 1.60218 \cdot 10^{-19}</math> Дж – электрон-вольт.
Минимальный период колебаний решетки в этих единицах – величина
порядка единицы. Функции, осуществляющей один шаг интегрирования
по времени, передается массив частиц и бокс. Имеются две
разновидности функций, которые можно использовать. Они определяют
вид граничных условий: условие отражения атомов от стенок бокса и
периодические граничные условия по двум направлениям в
пространстве. В экспериментах с пластинами использовалось второе
граничное условие, чтобы пластина фактически не имела свободных
краев.

Есть некоторая сложность, связанная с таким подходом. Дело в том,
что пространственный бокс имеет фиксированный размер. Это означает,
что при нагревании пластина переходит в напряженное состояние из-за
тепловых расширений. С другой стороны, такой подход позволяет задать
деформированную (растянутую) решетку в начальный момент времени.
Так как тепловые расширения малы, при нагревании напряжения в пластине
остаются напряжениями растягивающими. Эксперименты с атомами на поверхности
пластины показали, что деформации растяжения 3 процента практически
не влияют на поведение этих атомов.

Функция, осуществляющая один шаг интегрирования по времени, имеет
следующую структуру. В начале пространственный бокс разбивается на
кубические области. Размер ребра куба не меньше радиуса обрезания
<math>R + D</math>, определенного потенциалом взаимодействия. Это означает,
что каждый атом внутри одного куба может взаимодействовать только
с атомами внутри этого куба и атомами соседних кубов. Далее все
атомы сортируются по вышеописанным кубическим ячейкам. Чтобы
обеспечить условие периодичности, к граням бокса, на которых
поставлено это условие, добавляется по одному слою кубических
ячеек. В эти ячейки помещаются копии атомов из слоя кубических
ячеек, примыкающих к противоположной грани бокса. Затем следует
основной цикл – перебор всех частиц, находящихся в боксе, в
котором осуществляется приращение скоростей частиц за один шаг
интегрирования под действием сил. Для каждой частицы <math>(i)</math> определяется
набор частиц <math>(j_{n})</math>, с которыми она взаимодействует, исходя из
пространственного разбиения и радиуса взаимодействия. Затем,
внутри этого цикла вызывается функция, осуществляющая приращения
скоростей частиц <math>(i)</math> и <math>(j_{n})</math>, обусловленные слагаемыми
потенциала <math>V_{ij_{n}}</math>. Последняя функция производит вычисления
характерные для потенциала Терсоффа. Представление энергии в виде
<math>E = \sum_i E_i = \frac{1}{2} \sum_{i, j (\neq i)} V_{ij},</math>
является общим для всех потенциалов, поэтому функция,
осуществляющая шаг интегрирования по времени фактически не
привязана к способу взаимодействия. На завершающем этапе
рассматриваемая функция вычисляет приращения координат.

На этапе написания программа отлаживалась с использованием
[[Потенциал Морзе|потенциала Морзе]]. Потенциал Морзе не относится к многочастичным
потенциалам. Его вид отличается от выражения для <math>V_{ij_{n}}</math>
[[Потенциалы Терсоффа, Бреннера|потенциала Терсоффа]] тем, что
вместо коэффициента <math>b_{ij}</math>, вносящего многочастичность, используется
константа. Задаваемая равновесная конфигурация в виде решетки
структуры алмаза оказывалась неустойчивой, и атомы со временем
принимали структуру, близкую к плотноупакованной решетке.

== Вычисление температуры ==
Количественное выражение для температуры в
экспериментах связано с представлениями об идеальном газе, в
котором потенциальные взаимодействия отсутствуют. Равновесие
рассматриваемой атомной системы с идеальным газом предполагает
выравнивание средних кинетических энергий, приходящихся на один
атом. Поэтому температура в рассматриваемой системе считается
связанной с кинетической энергией так же, как она связана в
одноатомном идеальном газе. Коэффициент пропорциональности равен
<math>2 / (3k)</math>, где <math>k</math> – постоянная Больцмана.

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

== Нагрев кристалла ==
[[Файл:image091.png|thumb]]
На этапе отладки и тестирования программы, использующей
вышеописанный метод, проводились эксперименты, связанные в
основном с нагреванием и охлаждением кристаллической решетки.
Расчет проводился для пластины с поверхностью 1,1,1 толщиной 6
бислоев с периодическими граничными условиями на торцах. Подача
энергии в систему производилась равномерно по времени, т.е. полная
энергия системы линейно изменялась со  временем. При этом
отслеживалось изменение температуры.

Ниже приведены графики изменения кинетической (1) и полной (2)
энергии, значения на вертикальной шкале переведены в Кельвины.
Виден горизонтальный участок кривой при 1800 К, на котором,
видимо, происходит плавление.

==  Свободная поверхность 1,0,0 ==
[[Файл:image092.png|thumb|Рис.1. Вид кристалла Si. Идеальная решетка (перед нагреванием).]]
[[Файл:image094.png|thumb|Рис.2. Вид кристалла Si после нагревания.]]
Ниже приведены результаты компьютерного эксперимента по нагреванию
кремниевой пластины с поверхностью 1,0,0. Верхняя грань кристалла,
изображенная на рисунках, является свободной (как и нижняя, не
видная на рисунках). На боковых гранях поставлены граничные
условия периодичности. Пластина была подвергнута нагреванию,
выдержке и охлаждению в течение одинаковых интервалов времени,
равных 1000 единиц времени (см. "Описание компьютерной модели").
Нагревание проводилось с помощью сообщения случайных приращений
скоростей атомам на каждом шаге интегрирования, охлаждение – с
помощью введения диссипативных внешних сил, действующих на каждый
атом. На рисунках изображены состояния до и после нагревания до 700 К.

На первом рисунке видны ряды атомов поверхности пластины,
соответствующие идеальной решетке. На втором рисунке можно
наблюдать результат образования димерных рядов, что
свидетельствует о неустойчивости идеальной поверхности в данной
модели. Таким образом, проведенные расчеты качественно описывают
поведение поверхности 1,0,0, наблюдаемое в реальных экспериментах.

== Атом на поверхности 1,1,1 ==
Рассмотрим теперь результаты численного эксперимента с адатомом на
поверхности 1,1,1. Для уменьшения объема вычислений и для простоты
наблюдения  за адатомом пространственный бокс взят меньшего
размера, чем в эксперименте с пластиной с поверхностью 1,0,0.
Граничные условия те же – периодичность по двум направлениям в
плоскости пластины (эти два направления и задают плоскость
пластины, поскольку размеры кристалла примерно одинаковы по всем
трем направлениям, как видно из рисунков). Свободные поверхности
на рисунках изображены сверху и снизу. Как видно, толщина пластины
– 5 бислоев.

Атом помещается на некотором расстоянии от поверхности кристалла,
как показано на первом рисунке. В начальный момент времени
кристалл имеет температуру 0К. Атому сообщается скорость по
направлению к поверхности кристалла. Затем кристалл постепенно
нагревается. На рисунке 2 видно, что адатом задержался над атомом
поверхностного бислоя, заполняя объемную решетку следующего слоя
(т.е. в эпитаксиальном положении). При нагревании выше 300К адатом
переходит в другое положение, в котором взаимодействует с тремя
атомами поверхностного бислоя. Положение адатома над поверхностью
существенно ниже, на рис. 3, 4 стрелки указывают на адатом.
[[Файл:image100.png|thumb]]

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

[[Файл:vtsaplin_p1.png|thumb]]
[[Файл:vtsaplin_p2.png|thumb]]

Приведенные рисунки относятся к задаче на 1 процессоре со 100
поверхностными атомами. Эксперимент проводился на 50 процессорах,
т.е. общее число адатомов – 5000.

Ниже изображены распределения адатомов по их смещению
относительно начальных положений при температуре 1888 К. Серии
1, 2, 3 диаграммы соответствуют 300, 600 и 900 единицам времени.

[[Файл:vtsaplin_file0.png]]

Одна ячейка на диаграммах соответствует 1/10 межатомного
расстояния.

Первый максимум на рисунках соответствует адатомам, не перешедшим
в другое положение, второй максимум – перешедшим. Поскольку в
начальный момент времени равновесные положения адатомов
задавались приближенно, исходя из предположения, что все связи
между атомами равны между собой, то первый максимум распределения
находится не в нуле. Все адатомы сместились вдоль оси,
перпендикулярной поверхности кристалла. На диаграммах видно, что
со временем количество смещенных частиц возрастает.

Ниже изображены графики изменения относительного количества
смещенных адатомов <math>(n/N)</math> со временем при различных
температурах. Смещенными считаются адатомы, перешедшие минимум
распределения по смещению, который находится между первым и
вторым максимумами.

[[Файл:vtsaplin_dif.png]]

== Оценка частоты перехода адатомов в смежные равновесные положения ==
По значению количества смещенных адатомов определяется частота
перехода адатомов в смежные равновесные положения. Вероятность
перехода оценивается по следующей формуле:

<math>P = 1 - \exp(- t / \langle t \rangle),</math>

где <math>\langle t \rangle</math> – среднее время перехода.
Соответственно, среднее время перехода вычисляется по формуле:

<math>\langle t \rangle = \frac{t}{|\log (1 - n / N)|},</math>

где <math>t</math> – время эксперимента, <math>n</math> – число перешедших атомов, <math>N</math>
– полное число поверхностных атомов. Поскольку число перешедших
атомов мало по сравнению с полным, то приближенно можно написать:

<math>\langle t \rangle = \frac{t}{(n / N)}</math>

Или, для частоты перехода:

<math>\omega = \frac{1}{\langle t \rangle} = \frac{(n / N)}{t}</math>

Т.е., частота перехода характеризуется наклоном прямой
относительного количества смещенных адатомов со временем. Наклон
вычисляется по приведенным выше графикам. Причем первое значение
(при <math>t = 0</math>) отбрасывается, поскольку из-за начальных возмущений
наклон кривой вначале имеет большее значение, как видно на графиках.

График изменения частоты перехода от температуры имеет вид:

[[Файл:vtsaplin_Fraq.png]]

Т.к. все характеристики процессов с температурной активацией
пропорциональны экспонентам типа <math>\exp(-E/kT)</math>, где <math>E</math> –
энергия активации процесса, <math>k</math> – постоянная Больцмана, <math>T</math> –
абсолютная температура, то построим эту характеристику в
аррениусовых координатах, т.е. зависимость логарифма частоты
перехода, или среднего времени перехода, от <math>(1000/T)</math>. В таких
координатах должна получиться прямая, из наклона которой определим
энергию активации процесса.

[[Файл:vtsaplin_Arren.png]]

Поскольку точка плавления кремния находится в рассматриваемом
диапазоне температур, следует вычислить два значения энергии
активации: для кривых до и после плавления. Энергия активации до
плавления составила <math>E = 2.017 \cdot 10^{-19}</math> Дж <math>=
1.259 \,\mbox{eV}</math>. Энергия активации после
плавления составила <math>E = 1.775 \cdot 10^{-19}</math> Дж <math>= 1.108
\,\mbox{eV}</math>. Погрешность вычислений этих величин равна 2 и 6
процентов соответственно. Погрешность рассчитывается как
среднеквадратичное отклонение от теоретических зависимостей,
параметры которых вычисляются методом наименьших квадратов.

== Оценка частоты перехода адатома в смежные равновесные положения при наличии смежного адатома ==
Определим частоту перехода адатомов на поверхности при наличии
адатома в смежном положении. Поместим на поверхность кремния пары
адатомов, после чего проведем тот же компьютерный эксперимент,
что и для одиночных атомов. Положение адатомов указаны на
рисунке. Полное число адатомов <math>N = 10000</math>.

[[Файл:vtsaplin_p3.png|thumb]]

Приведем графики изменения относительного количества
смещенных адатомов <math>(n/N)</math> со временем.

[[Файл:vtsaplin_dif_2.png]]

Графики изменения частоты перехода от температуры для одиночных
адатомов и пар адатомов имеют вид:

[[Файл:vtsaplin_fraq2.png]]

В аррениусовых координатах:

[[Файл:vtsaplin_Arren2.png]]

Энергии активации и погрешности для одиночных адатомов и пар
адатомов приведены в таблице.

{| class="wikitable"
| 
| 1 адатом
| 2 адатома
|-
| Энергия активации до плавления, eV
| 1.259
| 1.093
|-
| Погрешность
| 2 %
| 3 %
|-
| Энергия активации после плавления, eV
| 1.108
| 1.031
|-
| Погрешность
| 7 %
| 15 %
|}

== Литература ==

1. J. Tersoff,
New Empirical Model for the Structural Properties of Silicon, Phys.
Rev. Lett., Vol. 56, Num. 6 (1986)

2. J. Tersoff,
New empirical approach for the structure and energy of covalent
system, Phys. Rev. B, Vol. 37, Num. 12 (15 April 1988-II)

3. Sakir Erkoc,
Empirical many-body potential energy functions used in computer
simulations of condensed matter properties, Physics Reports 278
(1997) 79-105