Документ взят из кэша поисковой машины. Адрес оригинального документа : http://genphys.phys.msu.su/rus/diploma/diploma2015/Drozdov_OV.pdf Дата изменения: Fri Feb 5 20:08:11 2016 Дата индексирования: Sun Apr 10 00:53:48 2016 Кодировка: Windows-1251

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ ?МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М.В.ЛОМОНОСОВА?

ФИЗИЧЕСКИЙ ФАКУЛЬТЕТ

Кафедра общей физики

Дипломная работа
?Теоретическое исследование формирования золотых наноконтактов при постоянной температуре?

Выполнил студент 605 группы: Дроздов Олег Вячеславович

Научный руководитель: к.ф.м.н. С. В. Колесников Допущен к защите Зав. кафедрой

МОСКВА 2015


Оглавление

Введение Глава 1. Теоретические и экспериментальные исследования металлических наноконтактов Постановка задачи Глава 2. Методика компьютерного моделирования наноструктур при постоянной температуре

3

4 8

9

2.1 2.2

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

9

Термостаты . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Масштабирующий термостат . . . . . . . . . . . . . . . . . . . . 11 2.2.2 Термостат Берендсена . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.3 Термостат Андерсена . . . . . . . . . . . . . . . . . . . . . . . . . 12

Глава 3. Результаты и обсуждение

13

3.1 3.2 3.3

Моделирование канонического ансамбля с использованием различных термостатов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Растяжение золотого наноконтакта вдоль направления [110] . . . . . 16 Растяжение золотого наноконтакта с горизонтальным смещением . . 16 3.3.1 Смещение вдоль направления [110] . . . . . . . . . . . . . . . . . 16 3.3.2 Смещение вдоль направления [001] . . . . . . . . . . . . . . . . . 19

3.4

Растяжение золотого наноконтакта с вращением . . . . . . . . . . . . 21
23 24 31

Результаты и выводы Приложение Литература

2


Введение

В современном мире трудно переоценить значимость электронных устройств для хранения, воспроизведения и передачи информации. Одним из основных направлений в физике конденсированного состояния вещества является теоретическое исследование возможности создания структур для производтсва более ?мких накопителей и увеличение скорости их работы. Использование металлических наноструктур открывает путь для развития подобных устройств. Существуют различные наноструктуры, претендующие на их будущее применение в электронике, такие как: нанокластеры, нанопровода, наноконтакты и тд. Данная работа посвящена изучению формирования наноконтактов. Для определения оптимальных условий для изготовления наноконтактов требуется детальное теоретическое исследование диффузионных процессов на атомном уровне. При исследовании образования контакта механическим способом, одним из важнейших с практической точки зрения является вопрос о влияния таких факторов, как скорость и направление растяжения. Экспериментально было установлено, что при растяжении вдоль определ?нных направлений контакт формируется занчительно более длинным и стабильным. Настоящая работа посвящена исследованию формирования золотых наноконтактов при постоянной температуре в широком диапазоне динамических условий растяжения. Определяются наиболее подходящие скорости и направления растяжения.

3


Глава 1 Теоретические и экспериментальные исследования металлических наноконтактов
Начало исследованиям атомных контактов было положено в работе [1], где впервые было описано получение золотых наноконтактов, имеющих проводомисть равную минимально возможному значению (кванту проводимости G0 =
2e2 h

). Они су-

ществовали ограниченное время при достаточно низкой температуре. Статья так же примечательна тем, что в ней впервые был предложен экспериментальный метод определени структуры атомных контактов и показана возможность формирования атомных наноконтактов длиннее одного атома. Существует несколько спобосов формирования наноконтактов, в их числе СТМ (скинирующая туннельная микроскопия), АСМ (атомно-силовая микроскопия), ПЭМ (просвечивающая электронная микроскопия) и метод механически управляемого разрыва. Методы АСМ и СТМ позволяют не только создавать, но и исследовать атомные контакты. С помощью СТМ уда?тся создавать нанопровода длиной до 10 атомов всего за несколько секунд. Одновременно с этим данный метод позволяет определять проводимость контакта, которая при растягивании уменьшается дискретно с шагом G0 и выходит на плато при формировании линейного контакта [2] (Рис.1). Однако, не всю информацию о формировании контакта можно получить экспериментальными методами. Некоторые вопросы может прояснить только теоретическое исследование, моделирование. Например, теоретически [3, 4] было установлено, что удлиннение нанопровода происходит за сч?т разрыва связей внутренних атомов, а в формировании самого контакта участвуют атомы, находящиеся на поверхности. Очень важным фактором при формировании наноконтакта является температура процесса. При низких температурах не активизируются диффузионные процессы, 4


G Рис. 1: Относительная проводимость G0 в зависимости от смещения (d) двух электродов при формировании золотого наноконтакта при температуре 4.2 K.

позволяя контакту дольше оставаться в растянутом состоянии, в то же время он становится более хрупким и дальнейшее воздействие может его разорвать. При повышении температуры контакт становится более тягучим, но поверхностная диффузия атомов может приводить к неустойчивости атомной цепочки. На примере меди было установлено, что контакт формируется с меньшим напряжением при температуре 200 К [5]. В современных теоретических работах производятся попытки более детально смоделировать эксперимент и найти причины несовпадения теории и эксперимента в классических работах. Например в статье [6] авторы смещают один из электродов в горизонтальном направлении (Рис. 3). Не менее важен вопрос о том, какие металлы позволяют формировать атомные контакты лучше. Чтобы это оценить измеряются различные физические величины. Например, в одной из статей предлагается оценивать отношения силы разрыва в бесконечной цепочке к силе разрыва кристалла [7] . В более поздней работе [8] в роли параметра для оценки вводится отношения энергии связи в цепочке к энергии связи в кристалле. Опираясь на эти критерии было установлено, что при механическом подходе к формированию контактов лучшим материалом является золото, после него идут иридий и платина, в то время как серебро и медь плохо подходят для создания атомных наноконтактов.

5


Рис. 2: . Зависимость механического напряжения (S) наноконтакта от растяжения (d) при различных температурах: (А) 100К, (Б) 200К, (В) 300К .

6


Рис. 3: Начальные конфигурации и направелния растяжения, синим цветом обозначены закрепл?нные атомы, ж?лтым свободные.

7


Постановка задачи

Основной целью данной работы является определение диапазонов скоростей и наравлений растяжения и вращения контакта, в которых он формируется более длинным и стабильным, в частности, были поставлены следующие задачи: 1. Разработать комплекс программ для реализации алгоритма молекулярной динамики, дополнить е? термостатами для моделирования систем при постоянной температуре. 2. Провести моделирование растяжения наноконтактов в широком диапазоне внешних условий и выявить зависимости между результатами моделирования контакта и внешними условиями растяжения.

8


Глава 2 Методика компьютерного моделирования наноструктур при постоянной температуре
2.1 Метод молекулярной динамики

Основным методом, использованным в данной работе, является метод молекулярной динамики [9], известный ещ? с середины прошлого века. Появления этого метода относят к пионерским работам Олдера и Вайнрайта [10], предложивших в 1957 году метод вычисления фазовых траекторий систем тв?рдых шаров. Классический метод молекулярной динамики исходит из сопоставления реальному молекулярному объекту математического образа системы материальных точек, движение которых описывается системой уравнений Ньютона:

d2 r i = Fi (1) dt2 Если рассматривать малый интервал времени (dt), то силу F можно считать постоянной. Величина интервала времени должна быть значительно меньше любого из периодов движения, свойственного системе. В данной работе dt = 1 Ч 10-14 c = 10 фс. Значение силы Fi (r1 , ...rn ) , зависящей только от координат, можно рассчитать, зная потенциальную энергию U (r1 , ...rn ) системы: m
i

Fi = -

U (r1 , ...rn ) ri

(2)

Потенциалы взаимодействия берутся в форме Клери-Росато [11]:

EC =
i

i ER + E

i B

(3)

9


i EB

=-
j

2

exp -2q

rij
r0

1/2

-1 rij


(4) (5)

i ER = j

A1

rij
r0

-1 +A0 exp -p

r

0

-1

Для расч?тов использовались параметры межатомного взаимодействия, представленные в работе [12]. Решение системы дифференциальных уравнений первого порякда может быть получено с помощью применения оператора эволюции:

(t) = e

iLt

(0)

(6)

iL =

(7)

где iL - оператор Лиувилля и Г многомерный вектор независимых переменных(координат и скоростей). Обычно, результат действия оператора эволюции на координаты не может быть посчитан аналитически. В таком случае используется разложение в ряд по малому параметру t = t/P , которое применяется P раз:
P

(t) =
i=1

(
s

eiL

s

t

)(0)

(8)

и вносит ошибку t /P
n

n-1

при разложения до n.

Используя гамильтонов формализм, оператор эволюции Лиувилля записывется следующим образом:
N N

iL =
i=1

v

i

ri

+
i=1

Fi (r) mi

v

i

(9)

Применяя формулу Троттера [13]:

eA

+B

= lim

N inf

eN e

A

B N

N

(10)

где A и B необязательно коммутирующий операторы, получим разложение оператора эволюции:

eiL
здесь iL1 =
N i=1 Fi (r) mi v
i

t

= eiL

t 12

eiL

2

t iL

e

t 12

(11)

vi ri В явном виде приближ?нный оператор эволюции, позволяющий получить коор10

и iL2 =

N i=1


динаты и скорости за время t записывается следующим образом:

r(t) = r(0) + tv(0) +

t2 F[r(0)] 2m

(12) (13)

t2 (F[r(0)] + F[r(t)]) 2m Таким образом мы приходим к скоростному методу Верле [14]. v(t) = v(0) +

2.2

Термостаты

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

kB T0 =
где

2 , s
tcalc

(14)

< K >=

1 tcal
c 0 N

K (t)dt, mi vi 2 (t) , 2

(15) (16)

K (t) =
i=1

s = 3N - 6 полное число степеней свободы системы. Для поддержания постоянной температуры модельной системы необходимо взаимодействие с внейшней средой, термостатом. При этом уравнения движения модифицируются определ?нным образом.
2.2.1 Масштабирующий термостат

Простейшим способом сохранения кинетической энергии системы является использование масштабирующего термостата [15]. Для этого скорости умножаются на фактор S = (
Tref 1/2 ) T

, где Tref требуемая температура. После чего система релакси-

руется некоторое число шагов, далее процедура повторяется.
2.2.2 Термостат Берендсена

Поместим теперь модельную систему в вязкую жидкость. Тогда уравнения Ньютона можно заменить на уравнения Ланжевена [16], т.е. рассматривать движение атомов как движение броуновских частиц:

dri, = vi, dt
11

(17)


mi
где

dvi =F dt

i,

- v

i,

+

Ri, (t) mi

(18) (19)


< Ri, (t) >= 0 < Ri, (t)Rj, (t + ) >= 2mi kB Tref ( )ij
основе его работы лежат следующие соотношения:

(20)

Следствием ланжевеновской динамики являеться термостат Берендсена [17]. В

dK 1 = lim t 0 t dt

i,

1 2 mi vi, (t + t) - 2

i,

1 2 mi vi, (t) 2

(21)

После элементарных преобразований и учитывая формулу для дисперсии Ri, , получаем:

dT dt

= 2 (Tr

ef

- T)

(22)

Что равносильно дополнительному слагаемому в уравнениях Ньютона:

dri, = vi, dt mi dvi =F dt
i,

(23) (24)

-

Tref - 1 vi, T

Параметр является характеристикой термостата, обратным временем релаксации.
2.2.3 Термостат Андерсена

Рассмотрим метод моделирования при постонной температуре, предложенный Андерсеном [18]. Связь системы с термостатирующей средой осуществляется за сч?т стохастических столкновений, которые испытывают случайные атомы. Между столкновениями система являается микроканоническим ансамблем и движется согласно обычным законам Ньютона. Стохастические столкновения гарантируют, что система не нарушает распределение Гиббса. Связь системы с термостатом характеризуется частотой случайных столкновений . Последовательные столкновения не коррелированы, следовательно распределения интервалов времени между ними является Пуассоновым [19, 20].

P (t, ) = e

- t

(25)

Где P (t, )dt вероятность того, что следующий удар произойд?т в интервале от t до t+dt. 12


Глава 3 Результаты
3.1 Моделирование канонического ансамбля с использованием различных термостатов
Прежде чем перейти непосредственно к исследованию факторов влияющих на формирование наноконтакта, необходимо определить, какой из вышеописанных термостатов больше подходит для решения данной задачи. Тестирование различных термостатов происходило на двух модельных системах. Первая представляла собой поверхность золота (001), которая релаксировалась при температуре 300 К. В рамках первой модели термостаты показали себя одинаково хорошо, система достаточно быстро ( 1000 шагов) достигала требуемой температуры и сохраняла е? постоянной. Было проверено, нарушают ли используемые термостаты распределения Максвелла по модулю скорости, и установлено, что лучше всего распределению Максвелла отвечает термостат Берендсена (Рис. 4). В качестве второй системы была выбрана модель наноконтакта, приближенная к экспериментальной. Из золота вырезалась цилиндрическая заготовка длинной 9 нм и диаметром 3.6 нм, ось симметрии которой совпадает с направлением [110]. В середине заготовки был произвед?н надрез, чтобы локализовать область формирования контакта. По два слоя сверху и снизу были закреплены. Термостатировались только верхние и нижние десять слоев для того, чтобы исключить влияние термостатов непосредственно на область образования наноконтакта. Моделирование происходило при постоянной скорости растяжения в vz = 2 Ч 10-1 м/с по оси z. В рамках тестовых расч?тов было установлено, что масштабирующий термостат и термостат Андерсена не подходят для моделирования таких систем. При их использовании возникает большая разность скоростей в приграничных слоях между термостатируемой и нетермостатируемой областью, вследствие чего именно там происходит разрыв заготовки. При использовании термостата Берендсена этого не происходит, так как в н?м 13


Рис. 4: Распределение количества частиц (N) по модулю скорости (V) при моделировании с использованием масшабириующего термостата (А), термостата Берендсена (Б) и термостата Андерсена (В). Сплошной линией обозначено теореитическое распределение Максвелла для данной температуры.

14

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

3.2

Растяжение

золотого

наноконтакта

вдоль

на-

правления [110]
На первом этапе исследования формирования золотых наноконтактов заготовка подвергалась растяжению в направлении [110], без смещения. Было установлено, что непосредственно перед началом образования одноатомного контакта происходит формирование структуры, практически сохраняющей кристаллическую фомрму, длинной свыше 1 нм и шириной в 2-3 атома в области разрыва (Рис. 5А). Эти результаты хорошо согласуется с экспериментальными данными [21]. При дальнейшем растяжении вышеописанная структура начинает утончаться, образуя наноконтакт толщиной в один атом и длинной до нескольких атомов, после чего происходит разрыв (Рис. 5Б). В момент разрыва контакта наблюдается скачок температуры в области его формирования (Рис. 5В) из-за увеличеня скоростей атомов в момент разрыва. В дальнейшем температура вновь выравнивается за сч?т термостатирования. Также в момент разрыва наноконтакта наблюдается скачок потенциальной энергии системы (Рис. 5Г) вследствие резкого уменьшения числа атомных связей в контакте.

15


Рис. 5: Структура наноконтакта:(А) в момент начала образования атомного контакта, (Б) непосредственно перед разрывом; (В) зависимость температуры (T) и (Г) потенциальной энергии (E) от удлиннения контакта (d). Потенциальная энергии наноконтакта отсчитывается от потенциальной энергии в момент разрыва.

3.3

Растяжение золотого наноконтакта со смещением

В реальном эксперименте очень сложно добиться смещения контакта в какомлибо конкретном кристаллографическом направлении. Это происходит в силу того, что проволока является макроскопической системой и малейшие дефекты могут повлечь появление деформаций сдвига и кручения, что приводит к неопредел?нности направления растяжения. Обычно, в теоретических работах этот факт не учитывается и исследуются только контакты растянутые вдоль направлений [100],[110] и [111]. В представленной работе была поставлена цель изучить этот вопрос более подробно, моделируя растяжения с дополнительным смещением. Модель, на которой проводились исследования растяжения со смещениями такая же, как в пункте 3.1. Дополнительные смещения производятся вдоль направлений [001] и [110].
3.3.1 Смещение вдоль направления [1 0]

1

При моделировании растяжения с дополнительным смещением менялся угол между направлением [110] и реальным направлением растяжения. В зависимости от этого угла строились изображения области образования контакта (Рис. 6). Угол 16


изменяется от 2 до 10 градусов с шагом 2 градуса. Исследование растяжения контакта со смещением в направлении [110] показало, что уч?т последнего сильно влияет на процесс формирования контакта. Уже при угле смещения в 2 градуса становятся заметны характерные особенности этого смещения. Во-первых, кристаллическая структура области образования контакта очень быстро разрушается, практически не образовывая устойчивого наноконтакта. При увеличении угла смещения не успевают образовываться пирамидальные структуры и основания контакта получаются более плоскими. Можно также отметить, что не происходит формирования переходной области растяжения контакта, как это происходило при растяжении без дополнительного смещения (Рис.5), так как заготовка раскалывается в самом нестабильном месте. Во-вторых, одноатомный контакт перед разрывом существует недолго что по-видимому связано с деформацией его структуры. Получившиеся в результате моделирования контакты похожи на экспериментально полученные при растяжении контакта вдоль направлений [100] и [111], которые являются невыгодными для образования наноконтактов. Одноатомный контакт вытягивается тем длиннее, чем больше угол дополнительного смещения, становясь при этом менее стабильным (Рис. 6).

17


18

Рис. 6: Структура нанопровода в момент начала образования атомного контакта (вверху) и непосредственно перед разрывом (внизу) в зависимости от угла смещения вдоль направления [110].


3.3.2

Смещение вдоль направления [001]

В случае дополнительно смещения вдоль направления [001] было установлено, что зависимость от угла смещения более плавная, чем в предыдущем случае. Процесс растяжения происходит без раскалывания заготовки. Формируется ярко выраженная переходная зона образования контакта, в которую со временем диффундируют атомы со всей поверхности заготовки. Эти атомы участвуют в формировании нанопроволок шириной в 2-3 атома, при этом их длина превышает 2 нм, а структура долго оста?тся похожа на структуру кристалла. При дальнейшем смещении деформация сдвига становится велика и атомные связи с противолежащей направлению смещения стороны начинают рваться. В результате образуется одноатомный наноконтакт длинной в несколько атомов, который сохраняется ограниченное время (Рис. 7). Длина одноатомного контакта зависит от начальной толщины области соединения. Последняя особенность присуща только этому направлению дополнительного смещения и не зависит от угла, проявляясь во всех модельных случаях. Более детальное моделирование момента разрыва показало, что такая структура наноконтакта более стабильна, чем в случае дополнительного смещения в направлении [110] и, в то же время, является менее протяж?нной, так как дополнительные атомы собираются в основании контакта, формируя пирамидальные структуры, не участвуя в его удлиннении. В случае продолжения растяжения оставшиеся связи рвутся, наноконтакт разрушается, оставляя пирамидальные структуры с вытянутыми в направлении контакта вершинами. Сравнивая полученные по растяжению с дополнительным смещением результатам можно придти к выводу, что растяжение происходит плавнее и образуются более длинные нанопровода в том случае, когда дополнительное смещение производится вдоль направления с наибольшим расстоянием между атомными слоями. При этом связи между атомами соседних сло?в слабее, что делает заготовку более тягучей и не происходит е? раскалывания во время растяжения. Также на образование контакта влияют процессы диффузии атомов по различным поверхностям кристалла. Исходя из этого, можно придти к заключению, что для различных целей использования наноконтактов можно использовать разные направления дополнительного смещения при растяжении: для стабильных широких нанопроводов лучше подходит направление [001], а для длинных одноатомных контактов направление [110].

19


20

Рис. 7: Структура нанопровода в момент начала образования атомного контакта (сверху) и непосредственно перед разрывом (снизу) в зависимости от угла смещения вдоль направления [001].


3.4

Растяжение золотого наноконтакта с вращением

При моделировании растяжения наноконтакта с дополнительным вращением изменялся угол между одинаковыми кристаллографическиими направлениями, перпендекулярными направлению [110], верхнего и нижнего электрода. В зависимости от этого угла строились изображения области образования контакта (Рис. 8). Угол изменялся в диапазоне от 2 до 10 градусов с шагом 2 градуса. В результате моделирования было установлено, что при вращении электрода в области контакта образуются структуры в виде двух пересекающихся пирамидок, характерные для растяжения контакта вдоль направления [111]. Так как это направление является невыгодным для формирования контактов, образование длинных атомных цепочек не происходит. Возможно, резкий обрыв наноконтактов обусловлен высоким механическим напряжением, возникающим из-за деформации кручения. Ещ? одной характерной особенностью этого варианта растяжения является смещение области формирования контакта от середины к закрепл?нному электроду. При угле закручивания выше 2 градусов на положения области разрыва больше влияет деформация кручения, нежели место предварительного утончения заготовки. Было установлено, что при дополнительном вращении область формирования контакта нестабильна и е? структура далека от кристаллической. Это приводит к тому, что между началом формирования контатка и его разрывом проходит достаточно мало времени. Сравнивая данные результаты и результаты по растяжению без дополнительных смещений и вращений можно сказать, что вращение вносит дестабилизирующий вклад в модельную систему, который приводит к ухудшению формирования атомных наноконтактов. Как следствие можно отметить тот факт, что для создания контактов экспериментально необходимо минимизировать число дефектов, которые могут привести к возникновению деформаций кручения.

21


22

Рис. 8: Структура нанопровода в момент начала образования атомного контакта (сверху) и непосредственно перед разрывом (снизу) в зависимости от угла поворота.


Результаты и выводы

1. Разработан программный комплекс реализующий метод молекулярной динамики при постоянной температуре. Данный комплекс позволяет решать широкий спектр задач, начиная от изучения фазовых переходов и заканчивая исследованием формирования различных наноструктур. 2. В качестве примера работоспособности программного комплекса было исследовано формирование золотого наноконтакта. При этом было установлено, что: (a) При дополнительном смещении одного из концов проволоки вдоль направления [001] вытагиваются стабильные нанопровода, имеющие толщину 2-3 атома и сохраняющие структуру близкую к кристаллической. В то время как при дополнительном смещении вдоль направления [110] образуются более длинные атомные наноконтакты. (b) При дополнительном закручивании контакта во время растяжения происходит разрушение кристаллической структуры в области контакта, вследствие чего атомные наноконтакты формируются очень плохо.

23


Приложение

В данном приложении представлена реализация метода молекулярной динамики на языке программирования Fortran90. Реализована возомжность закреплять верхние и нижние слои атомов, задавать скорость растяжения и угол вращения. Также реализованы следующие термостаты: Андерсена, Берендсена и масштабирующий. Блок отвечающий за последовательный вызов модулей программы:
call read_data call init_random_seed call maxwell do k = 1, nsteps call pos_vel if (tau.ne.0d0) then call thermo_berendsen else if (Rf.ne.0d0) then call thermo_andersen else if (iscal.ne.0d0) then call thermo_scale end if do i = natom - nstat_up + 1, natom call rotate end do t = 2 * ekin / Nf / kbt epot = er - ebound etot = epot + ekin end do

Модуль чтения входных данных
subroutine read_data open(1, file='Input.txt.')
24


open(2, file='Coords.txt') do i = 1, nparam read(1, *) parameter(i) end do do j = 1, natom read(2, *) coords(j), vels(j) end do end subroutine read_data

Модуль интегрирования уравнений движения
subroutine pos_vel do i = nstat_down + 1, natom - nstat_up x(i) = x(i) + vx(i) * h / 2 !x(i) = x(i) - sx * nint(x(i) / sx - 0.5) y(i) = y(i) + vy(i) * h / 2 !y(i) = y(i) - sy * nint(y(i) / sy - 0.5) z(i) = z(i) + vz(i) * h / 2 !z(i) = z(i) - sz * nint(z(i) / sz - 0.5) end do call force do i = nstat_down + 1, natom - nstat_up vx(i) = vx(i) + fx(i) * h / m vy(i) = vy(i) + fy(i) * h / m vz(i) = vz(i) + fz(i) * h / m x(i) = x(i) + vx(i) * h / 2 !x(i) = x(i) - sx * nint(x(i) / sx - 0.5) y(i) = y(i) + vy(i) * h / 2 !y(i) = y(i) - sy * nint(y(i) / sy - 0.5) z(i) = z(i) + vz(i) * h / 2 !z(i) = z(i) - sz * nint(z(i) / sz - 0.5) end do
25


do i = natom - nstat_up + 1, natom x(i) = x(i) + vx(i) * h y(i) = y(i) + vy(i) * h z(i) = z(i) + vz(i) * h end do return end subroutine pos_vel

Модуль составления списка парных атомов
subroutine pairs do i = 1, natom - 1 do j = i + 1, natom dxt = x(i) - x(j) dxt = dxt - sx * nint(dxt / sx) dyt = y(i) - y(j) dyt = dyt - sy * nint(dyt / sy) dzt = z(i) - z(j) dzt = dzt - sz * nint(dzt / sz) drt = sqrt(dxt**2 + dyt**2 + dzt**2) if ((drt.lt.Roff).and.(drt.gt.1e-7)) then np = np + 1 ip(np) = i jp(np) = j end if end do end do return end subroutine pairs

Модуль расч?та потенциальной энергии
subroutine energy do i = 1, np
26


w(i) = cut_off(dr(i), Ron, Roff) qexp1(i) = - (xi**2) * exp(- 2 * q * ((dr(i) / r0) - 1)) pexp1(i) = (a1 * (dr(i) / r0 - 1) + a0) * exp(- p * ((dr(i) / r0) - 1)) Eba(ip(i)) = Eba(ip(i)) - qexp1(i) * w(i) Eba(jp(i)) = Eba(jp(i)) - qexp1(i) * w(i) er = er + 2 * pexp1(i) * w(i) end do do i = 1, natom Eba(i) = sqrt(Eba(i)) ebound = ebound + Eba(i) end do return end subroutine energy

Модуль расч?ты сил
subroutine force call pairs call energy do i = 1, np dw(i) = cut_off_div(dr(i), Ron, Roff) qexp = -(q * xi**2) * exp(-2 * q * ((dr(i) / r0) - 1)) pexp = (a1 * p * (dr(i) / r0 - 1)+ a0 * p - a1)* exp(-p * ((dr(i) / r0) - 1)) f(i) = ((1 / Eba(ip(i)) + 1 / Eba(jp(i))) * (qexp + qexp1(i) * dw(i) / 2) + 2 * (pexp + pexp1(i) * dw(i))) / r0 / dr(i) f(i) = f(i) * 16 fx(ip(i))=fx(ip(i))+f(i)*dx(i) fx(jp(i))=fx(jp(i))-f(i)*dx(i) fy(ip(i))=fy(ip(i))+f(i)*dy(i)
27


fy(jp(i))=fy(jp(i))-f(i)*dy(i) fz(ip(i))=fz(ip(i))+f(i)*dz(i) fz(jp(i))=fz(jp(i))-f(i)*dz(i) end do return end subroutine force

Модуль, отвечающий за вращения атомов
subroutine rotate x_t = x - x0 y_t = y - y0 x = x0 + cos(d_thetta * pi / 180) * x_t + sin(d_thetta * pi / 180) * y_t y = y0 - sin(d_thetta * pi / 180) * x_t + cos(d_thetta * pi / 180) * y_t return end subroutine rotate

Модуль генерации распределения Максвелла
subroutine maxwell do i = nstat_down + 1, n - nstat_up, 2 call random_number (u1) call random_number (u2) x(i) = sigma * cos(2 * pi * u1) * sqrt(-2 * log(u2)) x(i + 1) = sigma * sin(2 * pi * u1) * sqrt(-2 * log(u2)) end do return end subroutine maxwell

Масштабирующий термостат
subroutine thermo_scale scale = sqrt(t_ref / t)
28


do i = nstat_down + 1, nthermo_down vx(i) = vx(i) * scale vy(i) = vy(i) * scale vz(i) = vz(i) * scale end do return end subroutine thermo_scale

Термостат Андерсена
subroutine thermo_andersen do i = nstat_down + 1, nthermo_down call random_number (R) if (R.lt.Rf) then call random_number (u1) call random_number (u2) vx(i) = sigma*cos(2*pi*u1)*sqrt(-2*log(u2)) * 1e10 call random_number (u1) call random_number (u2) vy(i) = sigma*cos(2*pi*u1)*sqrt(-2*log(u2)) * 1e10 call random_number (u1) call random_number (u2) vz(i) = sigma*cos(2*pi*u1)*sqrt(-2*log(u2)) * 1e10 end if end do return end subroutine thermo_andersen

Термостат Берендсена
subroutine thermo_berendsen lambda = 1 + h / (2 * tau) * (t_ref / t - 1) do i = nstat_down + 1, nthermo_down vx(i) = vx(i) * lambda
29


vy(i) = vy(i) * lambda vz(i) = vz(i) * lambda end do return end subroutine thermo_berendsen

30


Литература
[1] H. Ohnishi, Y. Kondo, K. Takayanagi, Quantized conductance through individual rows of suspended gold atoms, Nature 395, 780 (1998) [2] T. Kizuka, Atomic conguration and mechanical and electrical properties of stable gold wires of single-atom width, Phys. Rev. B 77, 155401 (2008) [3] F. Sato et al., Computer simulations of gold nanowire formation: the role of outlayer atoms, Appl. Phys. A 81, 1527 (2005) [4] A.L. Klavsyuk et al., Molecular dynamics study of Co-Au and Ag-Au bimetallic atomic chain formation, Eur. Phys. J. B 85, 331 (2012) [5] Liu, F. Wang, J. Zhao, L. Jiang, M. Kiguchib and K. Murakoshib, Theoretical investigation on the inuence of temperature and crystallographic orientation on the breaking behavior of copper nanowire, Chem. Phys. 11, 6514 (2009) [6] F. Tavazza, L. E. Levine and A. M. Chaka, Elongation and breaking mechanisms of gold nanowires under a wide range of tensile conditions, J. Appl. Phys. 106, 04352 (2009) [7] S.R. Bahn,K.W. Jacobsen, Chain Formation of Metal Atoms, Lett. 87, 266101 (2001)
Phys. Rev.

[8] А.Л. Клавсюк, С.В. Колесников, Е.М .Смелова, А.М. Салецкий, Моделирование процесса формирования металлических наноконтактов методом молекулярной динамики, ФТТ 53, 2356 (2011)

31


[9] Д.В. Хеерман, Методы компьютерного эксперимента в теоретической физике, М.: Наука, 1990. [10] B.J. Alder, T.E. Wainright, Phase Transition for a Hard Sphere System, J.Chem. Phys. 27, 1208-1209 (1957). [11] F. Cleri, V. Rosato, Tight-binding potentials for tranzition metals and alloys, Phys. Rev. B 48, 22 (1993). [12] N.N. Negulyaev, V.S. Stepanyuk, W. Hergert, P. Bruno, and J. Kirschner, Atomic-scale self-organization of Fe nanostripes on stepped Cu(111) surfaces: Molecular dynamics and kinetic Monte Carlo simulations, Phys. Rev. B 77, 085430 (2008). [13] H.F. Trotter, On the product of semi-groups of operators, math. Soc. 10, 545 (1959).
Proc. Amer.

[14] W.C. Swope, H.C. Andersen, P.H. Berens and K.R. Wilson, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, J. Chem. Phys. 76, 637 (1982) [15] W.G. Hoover, A.J.C. Ladd and B. Moran, High-Strain-Rate Plastic Flow Studied via Nonequilibrium Molecular Dynamics, Phys. Rev. Lett. 48, 26 (1982). [16] И.А. Квасников, Термодинамика и статистическая физкика. Том 3. Теория неравновесных систем. М.: Едиториал УРСС, 2003. [17] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A.DiNola and J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81, 3684 (1984). [18] H.C. Andersen, Molecular dynamics at constant pressure and/or temperature. J. Chem. Phys. 72, 2384 (1980).

32


[19] W. Weller, An introduction to Probability Theory and Its Application, Vol. 1, Wiley, New York, 1957. [20] W. Weller, An introduction to Probability Theory and Its Application, Vol. 2, Wiley, New York, 1966. [21] V. Rodrigues, T. Fuhrer, D. Ugarte, Signature of Atomic Structure in the Quantum Conductance of Gold Nanowires, Phys. Rev. Lett. 85, 4124 (2000).

33