Компьютерное моделирование нанотехнологий, наноматериалов и наноструктур : моделирование наносистем методами молекулярной динамики

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РФ 

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ  
ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ  
«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТЕХНОЛОГИЧЕСКИЙ УНИВЕРСИТЕТ «МИСиС» 

 

 
 
 

 

 

 

 
 

 

№ 2208 

Кафедра полупроводниковой электроники и физики полупроводников 

С.Ю. Юрчук 
 

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

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

Курс лекций 

Рекомендовано редакционно-издательским 
советом университета 

Москва  2013 

УДК 539.219.3(075.8) 
 
Ю83 

Р е ц е н з е н т  
д-р физ.-мат. наук, проф. Ф.И. Маняхин 

Юрчук, С.Ю. 
Ю83  
Компьютерное моделирование нанотехнологий, наноматериалов и наноструктур : моделирование наносистем методами 
молекулярной динамики : курс лекций / С.Ю. Юрчук. – М. : 
Изд. Дом МИСиС, 2013. – 47 с. 
ISBN 978-5-87623-663-0 

Цель курса лекций – ознакомить студентов с методами моделирования 
процессов выращивания наноструктур. Даны основные представления о методе Монте-Карло, основанном на выборе случайных процессов. Представлены приложения метода Монте-Карло для моделирования взаимодействия 
частиц. Даны основные представления моделирования выращивания нанопленок методом молекулярной динамики, с помощью которого можно численно интегрировать классические уравнения движения и проследить траекторию движения атомов и молекул в некотором конечном временном интервале, не превышающем нано- или микросекунду. Показаны примеры результатов моделирования роста нанослоев. 
Предназначен для студентов, обучающихся в магистратуре по направлению 210100 «Электроника и наноэлектроника». 
УДК 539.219.3(075.8) 

ISBN 978-5-87623-663-0 
© С.Ю. Юрчук, 2013 

ОГЛАВЛЕНИЕ 

1. Моделирование наносистем методами Монте-Карло.......................4 
1.1. Генерация случайных чисел на отрезке [a, b)  
в соответствии с заданной функцией распределения P(x) ..............5 
1.2. Интегрирование методом Монте-Карло.....................................7 
1.3. Приложения метода Монте-Карло к наносистемам, 
состоящим из нескольких частиц ......................................................9 
1.4. Применение метода Монте-Карло к неравновесным  
задачам................................................................................................10 
1.5. Уравнение Ланжевена................................................................11 
1.6. Взаимодействующие системы...................................................14 
Выводы...............................................................................................16 
2. Моделирование наносистем методами молекулярной динамики....17 
2.1. Принципы молекулярно-динамического моделирования 
наносистем .........................................................................................17 
2.2. Интегрирование уравнения движения Ньютона .....................21 
2.3. МД моделирование систем в контакте с тепловой  
ванной-термостатом..........................................................................25 
2.4. Плазмохимическое осаждение, моделирование роста 
покрытий в условиях бомбардировки высокоскоростными 
атомами...............................................................................................30 
2.5. Потенциал межчастичного взаимодействия............................32 
2.6. Алгоритм Бимона .......................................................................35 
Библиографический список...................................................................46 
 

1. МОДЕЛИРОВАНИЕ НАНОСИСТЕМ 
МЕТОДАМИ МОНТЕ-КАРЛО 

Моделирование методами Монте-Карло можно использовать в 
нанонауке для воспроизведения различных сложных физических явлений, включая расчеты свойств, предсказания фазовых переходов, 
самосборку структур, усредненных с помощью термической обработки, распределения зарядов и др. Имеется большое многообразие 
методов моделирования Монте-Карло, применяемых в зависимости 
от рассматриваемой наносистемы и желаемых результатов вычислений. Они включают классический, квантовый и объемный методы, 
но не ограничиваются ими. 
Основу моделирования Монте-Карло составляют случайные числа. Для успешного моделирования самая важная проблема – это генерация случайных чисел. 
Случайные числа, равномерно распределенные на отрезке [0,1), 
могут быть сгенерированны с помощью следующего алгоритма: 
1. Задается начальное значение чисел g, m. 
2. Рассчитывается случайное значение Sk = (gSk–1 + i)mod(m). 
3. Для k > 0 и 0 ≤ Sk < m случайное число rk = Sk/m, где k = 1, 2, …, ∞. 
S0 называется зерном. При использовании одного и того же зерна 
может воспроизводиться одинаковый набор чисел. Символы g, i и m 
обозначают генератор (или множитель), приращение и числовую характеристику. Функция «mod» на языках С++ и Delphi возвращает 
остаток (на языке Visual Basic функция RND). Предпочтительно, чтобы S0 было большим нечетным числом. Число m целесообразно выбирать равным 2Е, где Е > 35. Большое число g уменьшает корреляции 
внутри ряда. Число i должно быть простым по отношению к числу m, 
g – 1 многократным по отношению к любому элементарному числу m. 
В этом случае период последовательности будет равен 2Е. В языках 
высокого уровня имеются встроенные функции генерации случайных чисел (Visual Basic и Delphi – «Random», для инициации генератора случайных чисел имеется функция «Randomize», которая требуется, чтобы исключить повторяемость набора случайных чисел). 

1.1. Генерация случайных чисел на отрезке [a, b) 
в соответствии с заданной функцией 
распределения P(x) 

Выборка по значимости, также известная как «смещенная выборка», может использоваться для генерации случайных чисел согласно заданной функции распределения P(x). Это возможно потому, 
что эффективность методов выборки по значимости может быть заметно увеличена соответствующим отбором точек выборки. Описываемая методика называется «выборкой по значимости». 
Назначение выборки по значимости состоит в том, чтобы заострять внимание на точках выборки, где функция распределения значительна, и избегать выборки там, где значение функции пренебрежимо мало. Необходимо внести поправку на это смещение точек выборки, взвешивая определенным образом значения выборки. Общий 
подход к выборке по значимости состоит из двух этапов. 
1. На начальном этапе проводится предварительная оценка функции P(x). Интересуемая область делится на довольно крупные ячейки, значение функции P(x) определяется по центру каждой ячейки, и 
это значение и кумулятивная сумма записываются. 
2. Следующий этап состоит в отборе стольких точек выборки, 
сколько желательно. Но применять методику отбора надо таким образом, чтобы число точек выборки в каждой ячейке было пропорционально значению функции P(x) в центре этой ячейки, определенному на начальном этапе. При этом каждой выборке задается 
вес, обратно пропорциональный значению функции в центре ячейки. 
Существует несколько способов генерации случайных чисел при 
помощи выборки по значимости. 
1. В случае заданного отрезка [a, b) для получения требуемого 
числа случайных точек N сначала разделим интервал на 
N  частей 
(это число выбирается произвольно). Далее создадим однородный 
набор точек для каждой части, при этом число генерированных точек 
будет равно 

 
( )
(
)
i

i

NP x

P x
∑
, 
(1.1) 

где xi – точка в середине каждой части. Очевидно, что генерированные точки распределяются в соответствии с Р, и их число почти N. 

Однако следует заметить, что положения точек сильно скоррелированы и перед использованием требуется их переставить. 
2. Однородные случайные числа в интервале [0, 1) отрезка [a, b) 
можно отобразить аналитически, если функция распределения Р(х) 
достаточно проста. Если функция распределения вероятностей r в 
интервале [0, 1) однородна, то имеем  

 
( )
( )
(
)
(
)
,
;
,
1.

x

a
dr
P x dx
r
P t dt
g x a
g a b
=
⇒
=
=
=
∫
 
(1.2) 

Это соотношение необходимо преобразовать, чтобы получить выражение для х на отрезке [a, b) в соответствии с Р равномерно. 
3. Третий способ состоит в использовании известного алгоритма 
Метрополиса для произвольно усложненной функции Р. В алгоритме 
можно начать с произвольной точки в интервале [a, b), например 
х0 ≈ (a + b)/2, затем добавить к х0 случайное число, согласно выражению xi = x0 + d(2r – 1), где d – шаг; r – случайное число в [0, 1). Если 
(
)
(
)
0
1
i
P x
P x
> , то шаг принимается; в противном случае это отношение сравнивается со случайным числом в интервале [0, 1): если 
данное отношение больше, тогда шаг снова принимается, иначе он 
отвергается, а предыдущая точка х0 сохраняется, и начинается другое 
испытание из этой точки. Имеется также вероятность выхода генерируемых точек из области значений [a, b); в этом случае эти точки отбрасываются. По сравнению с первым способом этот метод намного 
менее эффективен, потому что дает много выброшенных точек. Эффективность данного метода очень зависит от выбора размера шага 
d. При маленьком d большинство испытаний будет принято, но для 
хорошей выборки потребуется много испытаний. Если же d большой, 
то можно произвести выборку в [a, b) быстрее, но тогда может быть 
много выброшенных точек. 
4. Более эффективный способ, предложенный фон Нейманом в 
1951 г., основан на выборе точек функции f {для х на отрезке [a, b); 
f(x) > P(x)} с использованием представленного выше второго способа 
и генерации точек, распределенных согласно f (рис. 1.1). Наиболее 
просто было бы принять f постоянной, но большей, чем максимальное Р. Когда такая точка генерируется на отрезке [a, b), случайное 
число r в интервале [0, 1) сравнивается с отношением Р(х)/f(x). Если r 
меньше, чем это отношение, то число х принимается, в противном 
случае оно отвергается, и следует повторить этот процесс с другой 

точкой х, распределенной согласно функции f. Для постоянной f вероятность принятия х, очевидно, пропорциональна Р(х). 
Схема принятия гипотезы в алгоритме фон Неймана для выборки 
значений двух случайных чисел, распределенных в интервале [0, 1) 
(см. рис. 1.1), следующая: первое число, преобразованное по методу 
обратных функций с распределением f, нанесено по координате x, 
второе, преобразованное с равномерным распределением, по координате y между 0 и f(x). Если эти числа попадают в область значений 
под функцией P(x), то точка принимается. 

 

Рис. 1.1. Схема принятия гипотезы в алгоритме фон Неймана  
для выборки значений двух случайных чисел,  
распределенных в интервале [0, 1) 

1.2. Интегрирование методом Монте-Карло 

Стандартные методы численного интегрирования – методы прямоугольников, трапеций, Симпсона, квадратуры Гаусса. Эти методы, 
даже если они очень точны для одномерных (обычных) интегралов, 
быстро обращаются в предельно допустимые при расчетах многомерных (многократных) интегралов, так как необходимое число оценок функции растет весьма быстро, как 
,
d
N
 где N – число оценок 
функции, требуемое для вычисления одномерного интеграла, а d – 
размерность многомерного интеграла. В этой ситуации интегрирование методом Монте-Карло, использующее случайные числа, может 
стать очень удобным. 

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

 
( )
(
)

1

.

b
N

i

i
a

b
a
I
f x dx
f x
N
=

−
=
≈
∑
∫
 
(1.3) 

Точки хi случайно и равномерно распределены на отрезке [a, b). 
Так как число точек N растет, центральная предельная теорема 
(ЦПТ) гарантирует конвергенцию суммы к реальному значению интеграла, причем ошибка имеет порядок 1
.
N  
Вышеупомянутый метод интегрирования, хотя и сходящийся, повидимому, не очень эффективен, если на данном интервале функция 
имеет большие осцилляции. Рассмотрим, например, функцию 
( )

10x
f x
e
=
 в интервале [0, 1). Разумеется, большое число точек и 
оценки функции будут «бесполезными» вблизи х = 0, где значение 
функции намного меньше (более чем на 4 порядка), чем для точек 
вблизи х = 1. Выборку чисел на интервале лучше производить неравномерно, чтобы накопить больше точек вблизи х = 1 и меньше точек 
в области, близкой к х = 0. Это можно осуществлять выборкой по 
значимости. Подберем неоднородную функцию распределения Р(х), 
согласно которой будут выбираться точки хi. Кроме того, можно записать интеграл I в следующей форме: 

 
( )
( )
( )
( )
(
)
(
).

b
b
N

i
i
i P
a
a

b
a
I
f x dx
P x
f x
P x
dx
f x
P x
N
∈

−
=
=
⎡
⎤
≈
⎣
⎦
∑
∫
∫
 (1.4) 

Заметим, что теперь интегрирование функции f(x)/P(x) проводится 
в точках, распределенных согласно Р(x). Функция Р(x), конечно, аналитически известна и нормализована к единице. Флуктуация будет 
меньше в подынтегральном выражении I/P(x), если функция распределения Р(x), которую можно выбрать по желанию, будет похожа на 
исходную функцию f(x). 
Основная идея выборки по значимости – по крайней мере в том, 
чтобы подобрать функцию Р(x), похожую на f(x), поэтому функция 
I/P(x) почти не флуктурирует, и разброс результатов вычислений незначителен. Другими словами, значение интеграла определяется точно даже при относительно небольшом количестве точек. 

Формула интегрирования непосредственно обобщается в кратные 
интегралы 

 
(
)
(
)

1

, ,
,
,
,
,
,

N
G

i
i
i
i
G

v
f x y
z dv
f x y
z
N
=
=
∑
∫∫ ∫
…
…
…
 
(1.5) 

где 
G
v  – объем области интегрирования. Например, для двукратного 
интеграла с прямоугольной области интегрирования имеем 

 
(
)
(
)(
)
(
)

1

,
,
.

b d
N

i
a c

b
a
d
c
f x y dxdy
f x y
N
=

−
−
=
∑
∫∫
 
(1.6) 

1.3. Приложения метода Монте-Карло 
к наносистемам, состоящим  
из нескольких частиц 

Подход к интегрированию вариационным методом Монте-Карло 
используется для оценки многомерных интегралов. Это расчет основного энергетического состояния многочастичных систем Е0: 

 

*

0
*
H
dR

E

dR

Ψ
Ψ
=
Ψ Ψ
∫
∫
, 
(1.7) 

как, например, нескольких электронов в атоме, молекуле или квантовой точке. 
В этом уравнении Н – гамильтониан и  

 
(
)
1
2
0
,
,
, N
i
i

i
r r
r
c

∞

=
Ψ
=
Ψ
∑
…
 
(1.8) 

– пробная функция, где коэффициенты ci ряда нормализованы 

 

2

0

1.
i
i
c

∞

=
=
∑
 
(1.9) 

Пробная волновая функция содержит ряд вариационных параметров. Интеграл (1.7) можно переписать через плотность распределения вероятностей волновой функции 

2
Ψ . Это выражение составляет 

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

3
2
E
NkT
f
=
+
*1
(
)
(
)
1
1
1
1
3
,
2

N
r
r
N
N
N
NkT
dr
dr
r
r
dr
dr e−βϕ
=
+
ϕ
∫
∫
…
…
…
…
 (1.10) 

где N – число молекул; k – постоянная Больцмана; Т – температура. 

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

2
Ψ  в квантовом случае и 
(
)
1
N
r
r
e−βϕ
…
 в классическом случае. 
Затем вычисляется среднее значение «локальной энергии» 

 
( )
( )
( )
(
)
1
2
;
,
,
,
,
N
E R
H
R
R
R
r r
r
=
Ψ
Ψ
=
…
 
(1.11) 

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

1.4. Применение метода Монте-Карло 
к неравновесным задачам 

Метод Монте-Карло может использоваться для моделирования 
временных изменений состояния стохастической системы во времени, когда известны скорости переходов. Здесь скорости переходов не 
могут быть произвольными. Их выбор описывает реальный неравно
––––––––– 

1* Угловые скобки 
 обозначают оператор математического ожидания. 

весный физический процесс перехода, обусловливающий динамику 
системы, и должен быть представлен заранее. 
В физических процессах матрица перехода не формируется произвольно, а ее можно каким-нибудь образом вывести или аппроксимировать. Например, рассмотрим процесс диффузии атома на поверхности. Так как этот атом может перемещаться налево, направо, 
вверх или вниз (простая модель с вероятностью перехода, деленной 
на kT, зависящей от энергетических барьеров, которые атому требуется преодолеть), то можно определить четыре матрицы переходов 
для четырех процессов. Эти элементы представляют собой четыре 
фактора Больцмана, которые известны, если известны четыре потенциальных барьера. Именно это имеется в виду, когда говорится о 
«матрице переходов реального процесса, известного заранее». 
Типичные проблемы, решаемые с помощью этого метода, – проблемы вычисления транспортных свойств, например, расчет переноса носителей в электронных устройствах, касающийся решения 
уравнения переноса Больцмана, переноса энергии, массы, импульсов, 
переносов агрегатирования и роста. 

1.5. Уравнение Ланжевена 

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

 
( )

2

2
,
d y
dx
f t
dt
dx
μ
μ
= −
+
τ
 
(1.12) 

где μ – коэффициент трения, определяемый из закона Стокса;  
τ – нормировочный коэффициент. 

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