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

1.2. КВАНТОВОЕ МОЛЕКУЛЯРНОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ДИФФУЗИИ ПРИ СТРУКТУРНОМ ФАЗОВОМ ПЕРЕХОДЕ В АМОРФНОМ УГЛЕРОДЕ

Попов Александр Михайлович, профессор, профессор, д-р физ.-мат. наук, факультет вычислительной математики и кибернетики, Московский Государственный Университет им. М.В. Ломоносова, Москва, действительный член РАЕН по отделению «Прикладная математика и управление», член Американского Физического общества, popov@cs.msu.su

Шумкин Георгий Николаевич, разработчик ПО, канд. физ.-мат. наук, IBM Systems and Technology Group, Научно-технический центр IBM в России, Москва, Россия; georgiy-sh@yandex.ru

Никишин Николай Глебович, аспирант, факультет вычислительной математики и кибернетики, Московский Государственный Университет им. М. В. Ломоносова, Москва; mail@nikishinng.ru

Аннотация: Целью работы является моделирование, направленное на объяснение процессов, происходящих при создании памяти, основанной на фазовых переходах. В работе проводится квантовое моделирование «из первых принципов» процесса диффузии атомов в аморфном углероде при структурном фазовом переходе. Показано, что термические эффекты приводят к образованию графитовой послойной структуры в молекулярной системе. В окрестности термически обусловленного фазового перехода процесс диффузии является анизотропным и обусловлен взаимодействием ионов друг с другом с образованием электронных ковалентных связей. Временной процесс увеличения среднеквадратичного отклонения отличается от Эйнштейновской зависимости из-за взаимодействия атомов и образования связей. Такая структура соответствует повышению электрической проводимости по сравнению с аморфной структурой. Плотность носителей заряда увеличивается в окрестности графитовых поверхностей. Полученные характеристики диффузии позволяют объяснить механизм фазового перехода в экспериментах, в которых была показана возможность создания памяти на фазовых переходах на основе аморфного углерода [1].

Вычисления проводились на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ.

Ключевые слова: многомасштабные квантово-механические коды молекулярной динамики, фазовый переход в аморфном углероде, память на фазовых переходах, нанотехнологии, суперкомпьютер IBM Blue Gene/P

1.2. QUANTUM MOLECULAR MODELLING OF DIFFUSION PROCESS DURING PHASE CHANGE IN AMORPHOUS CARBON

Popov Alexander M., the Doctor of Science in Physics and Mathematics, professor of the Lomonosov Moscow State University, the full member of the Russian Academy of Natural Sciences on department "Applied mathematics and management", the member of the American Physical society, po-pov@cs.msu.su

Shumkin Georgiy N., PhD in Physics and Mathematics, Software Developer, IBM Systems and Technology Group, IBM Science and Technology Center in Russia, Moscow, georgiy-sh@yandex.ru

Nikishin Nikolai G., postgraduate, Computational Mathematics and Cybernetics Department, Lomono-sov Moscow State University, Moscow, mail@nikishinng.ru

Abstract: This paper purpose simulation performing to explain processes taking place in phase-change memory creation. Ab initio quantum simulation of atomic diffusion in amorphous carbon during phase change process is presented. It is shown that thermal effects lead to creation of a graphitic layers structure in a molecular system. Calculations show that in conditions of thermally induced phase change the diffusion process is anisotropic and forced by the fusion of atoms with different bonds and creation of covalent bonds. Time depended process of increasing mean square displacement differs from Einstein relation due to atoms interaction and formation of covalent bonds. Such a structure corresponds to increasing of electric conductivity in comparison with amorphous structure. The density of charge carriers is increased in the area of graphitic surfaces. Obtained dependencies help to understand the mechanism of phase change in experiments directed to the making of phase change memory on the base of amorphous carbon [1].

IBM Blue Gene/P supercomputer installed at the Faculty of Computational Mathematics and Cybernetics of the Moscow State University is used for calculations

Index terms: multiscale quantum molecular dynamics codes, phase transition in amorphous carbon, phase change memory, nanotechnology, supercomputer IBM BlueGene/P

1. Введение.

В последнее время много внимания уделяется использованию различных форм углерода, таких как нанотрубки и графен при создании энергонезависимой памяти, основанной на фазовых переходах (PCM) [1], [2]. По сравнению с этими формами углерода, аморфный углерод (a-C) в особенности интересен, так как легко позволяет включить себя в клеточную структуру наиболее удобную для технологии PCM[3]. Изменение атомной структуры материала, приводящее к переключению проводимости под действием внешних полей и температуры, лежит в основе создания памяти на фазовых переходах [4], [5]. В научно-исследовательской лаборатории IBM в Цюрихе была проведена серия экспериментов [1], в которых была показана возможность создания такой памяти на основе a-C.

Углерод существует в различных формах. Наиболее

2

часто встречающимися являются sp -доминирующая графитовая форма с низким сопротивлением и sp3 -доминирующая алмазная форма с высоким сопротивлением. При напылении пленки a-C имеет опре-

23

деленное отношение sp и sp атомов углерода. Изменение этого соотношения может существенно влиять на проводящие свойства углерода. Даже при

3 2 2

определенном отношении sp /sp кластеризация sp углерода может приводить к большому вкладу в проводимость [6], [7]. Элементный вклад углерода делает приборы, основанные на фазовых переходах, масштабируемыми на очень малые размеры. Изучение проблем диффузии атомов в a-C может помочь в решении проблемы уменьшения дрейфа сопротивления и повышения надёжности. Также изменение проводимости углерода по отношению к множеству стимулов делает материал привлекательным для использования в приложениях по созданию PCM.

В предыдущих исследованиях были изучены термические эффекты, приводящие к смене проводимости [8]. Также было показано, что сопротивление a-C может меняться при изменении тока, протекающего через материал [9].

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

В данной работе мы представляем численное моделирование атомных диффузионных процессов в a-C и их влияние на образование структур и смену сопротивления. Мы используем квантово-механические численные коды для исследования механизмов появ-

ления структуры и смены сопротивления проводящей наноточки.

Численному моделированию структур аморфного углерода посвящен целый ряд работ. В работах [10], [11] проводилось моделирование на основе классической молекулярной динамики с эмпирическими потенциалами для исследования различных структур a-C. Основным недостатком такого моделирования является многозначность появляющихся структур, что обусловлено неопределенностью в задании эмпирического потенциала взаимодействия между частицами. Квантово-механическое моделирование применялось для изучения структур тетраэдрального, алма-зоподобного аморфного углерода. Совместная модель связанных уравнений классической молекулярной динамики и транспортных уравнений сплошной среды предложена в [11] для описания памяти на фазовых переходах в аморфном углероде. Однако требуется дальнейшее исследование вопросов, связанных с согласованием поведения молекулярной системы на основе квантовой молекулярной динамики с пространственным распределением температуры, получаемой из уравнений сплошной среды. В работах [12] и [13] проведено компьютерное квантово-механическое моделирование процессов возникновения анизотропной атомной структуры в a-C. В них показано, что появление графитовой структуры из аморфного состояния обусловлено термическими эффектами. Тем не менее, механизм образования структуры не понят до конца. Повышение проводимости системы, связанное с изменением молекулярной структуры в эксперименте, должно приводить к существенному увеличению джоулева нагрева проводящей точки и дальнейшему повышению температуры [14].

В данной работе для моделирования атомной структуры используется квантово-механическая молекулярная динамика Кара-Парринелло, представленная в численном коде CPMD [15], [16], основанная на теории функционала плотности [17]. Рассматривается молекулярная система, представляющая аморфную структуру а-С. Основное внимание уделяется возникновению анизотропных параметров, характеризующих структуру молекулярной системы.

Расчеты в рамках квантовой молекулярной динамики требуют большой вычислительной мощности. Вычисления проводились на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ.

2. Квантовая модель молекулярной динамики и детали ее численной реализации.

При построении моделей молекулярной динамики из «первых принципов» исходят из уравнения Шре-дингера и кулоновского взаимодействия заряженных частиц. Мы используем математическую модель молекулярной динамики Кара-Парринелло [15], [16], [17] на основе приближения Борна-Оппенгеймера. В приближении Борна-Оппенгеймера считается, что движение массивных ядер (ионов) может быть описано в классическом приближении. Силы, действующие на ядра со стороны электронов, рассчитываются на основе решения уравнения Шредингера для электронов в основном состоянии. Мы будем использовать кванто-во-механическую формулировку, основанную на уравнениях функционала плотности (<^Т» [17]).

Уравнения движения ядер в усредненном поле электронов записываются в виде:

M

Z Z

R/(i) = -?VR, IRJ+jp(r)VR, IrZ-ri

Z

Dr

(1)

Для поиска новых орбиталей требуется решать уравнение Шредингера при каждом новом положении ядер:

-1À + U ff (р,r) \у (r) = е у (r) ■ ^ 2 eff in nTn

(2)

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

N

e 2

P(r)= Е (r) I2

(3)

п = 1

где величина Ne равна количеству валентных электронов в системе. Эффективный потенциал представляется в форме:

г, ггг, , (4)

и (р,г) = ин (р,г)+и (р,Ур,г)-у—-I— + У 1 ш н хс I IКI - г I I < у I КI - К ^ I

Где ин - потенциал Хартри, а ихс - обменно-корреляционный потенциал в форме РВЕ (Perdew-Burke-Ernzerhof) [15].

Электронные орбитали фп раскладываются по плоским волнам в приближении нулевой зоны Бриллю-эна:

(5)

Vn(r) = ЕCn,k exp(/k ■r)'

где О - объем ячейки, к - вектор обратной решетки, и сп,к - Фурье-коэффициенты орбитали п.

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

Приведем ряд деталей численной реализации. Она важна для понимания точности и возможностей численного моделирования. Для расчетов фазового перехода при помощи кода CPMD [15] использовался суперкомпьютер IBM Blue Gene/P, установленный на факультете ВМК МГУ. Самой затратной по времени вычислений частью является параллельный расчет электронной структуры системы. Используемый алгоритм обладает вычислительной сложностью, масштабируемой как O(Ne3).

Процесс решения задачи включает в себя нескольких этапов. Сначала по коэффициентам cn,k при помощи трехмерного быстрого преобразования Фурье вычисляются элементы матрицы гамильтониана H. Затем решается нелинейная задача на собственные значения, и проводится пересчет плотности и эффективного потенциала в реальном пространстве с использованием трехмерного быстрого преобразования Фурье. Для поиска собственных значений на каждой итерации в коде CPMD реализован параллельный алгоритм Ланцоша с использованием технологии MPI. В данной реализации используется быстрая раздельная сеть коллективных межпроцессорных сообщений суперкомпьютера IBM Blue Gene/P.

Размер суперъячейки, в которой располагались атомы углерода, составлял 14,268x10,701x10,701 Â. Начальная конфигурация включала 192 атома: 44 атома sp2 с тремя связями, представлявшие графитовые области; 78 атомов sp3 с четырьмя связями, соответствовавшие алмазным областям; 70 атомов sp с двумя связями, представлявшие линейные области. Соответствующее процентное содержание числа атомов в начальной конфигурации составляло 23%, 41% и 36%. Плотность системы была равна n = 2, 34 г/см3. При этом использовалось разложение по 384 электронным орбиталям, 54554 плоских волн для каждой волновой функции и 436234 плоских волн для плотности. Временной шаг был равен 0.125 фс. Пространственная сетка состояла из 144x108x108 точек.

3. Моделирование процесса диффузии атомов.

Ряд расчетов молекулярной динамики при разных температурах был проведен в работах [12], [13], [14]. Во всех расчетах при поднятии температуры алмазные области замещались графитовыми. Сравним установившиеся атомные структуры в процессе молекулярной динамики.

При температурах порядка T = 3000 K формируется графитоподобная послойная структура, обладающая меньшим сопротивлением, чем начальная конфигурация.

Эти изменения, обусловленные температурным эффектом, рассматриваются как фазовый переход в наноструктуре материала, приводящий к пороговой смене сопротивления. В ходе нагрева обнаруживается порог по температуре, выше которого возникают высокопроводящие кластеры графита.

Диффузия в аморфном углероде имеет сложный характер. Фазовый переход происходит при температуре Т = 2000 К. Атомы, обладающие различным числом связей - 5р, 5р2, sp3 - диффундируют по-разному. Для каждого сорта важно свободное пространство, в котором атомы могут разогнаться. Когда атом приближается на определенное расстояние к другому атому или атомной конфигурации, происходит их слияние в новую структуру. В зависимости от мобильности атома он может остановиться в своем продвижении и колебаться относительно нового положения. Мобильный атом с большой кинетической энергией может пройти сквозь всю элементарную ячейку.

На рис.1 (а, б, в, г) представлена диффузия в атомной конфигурации в виде четырех структур, полученных в четыре момента времени в процессе молекулярной динамики при Т = 3000 К. Разными цветами

2 3

отмечены атомы: sp красным, sp - серым, sp - синим. Хорошо видно, что все цвета кроме серого исчезают. В конце конфигурация представлена только графитовыми атомами. Таким образом, помимо смещения ионов происходит преобразование структур при изменении числа связей атомов. Возникновение высокопроводящих атомных конфигураций происходит в результате замещении алмазных областей графитовыми.

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

Структура, полученная в конце расчета молекулярной динамики, в среднем отражает параметры структуры графита. В работе [14] введен параметр структуры, характеризующий образование порядка из аморфной молекулярной структуры. Рассмотрена пространственная плотность ионов в процессе молекулярной динамики. Пространственная плотность ионов вычислялась с использованием техники гистограмм в процессе эволюции системы.

На основе статистики вычисленных траекторий, используя положения атомов при определенной температуре, рассчитывается структурный фактор S(t).

Изучение статистического процесса диффузии атомов проводится путем вычисления среднеквадратичного отклонения для каждой компоненты «а» (обо-

23

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

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

r Ца=N i(- ч о)

(6)

Рис. 1. Эволюция структуры молекулярной системы при температуре Т = 3000 К (а) аморфная конфигурация в начале процесса нагрева, (б), (в) атомная конфигурация в середине процесса нагрева (г) графитоподобная конфигурация в конце процесса. Синие атомы соответствуют атомам с $р3-2

связями, серые - атомам с $р -связями, красные - атомам с $р-связями.

Основываясь на этих данных можно оценить коэффициенты диффузии используя соотношение Эйнштейна для самодиффузии:

r.(t) - г.( о)2 } = 6Dt + C

(7)

где С - это константа и D - коэффициент самодиффузии.

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

Структурный параметр можно оценить следующим образом:

S (t) =1 (r 2(t)

где величина в скобках вычисляется как статистическое среднее по соответствующей группе атомов «а».

(7)

Диффузия атомов в аморфном углероде анизотропная, в особенности, если мы в начальных условиях в качестве затравки размещаем участки графитовых плоскостей. Характеристики процесса диффузии при Т = 3000 К в различных направлениях показаны на рис.2: (а) среднеквадратичное отклонение, определяемое по формуле (6), и выраженное в ангстремах в квадрате, и (б) отношение среднеквадратичного отклонения ко времени, выраженное в ангстремах в квадрате ко времени, выраженном в пикосекундах. Данное отношение должно сходиться по времени и, согласно соотношению Эйнштейна (если оно спра-

а

2

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

Рис. 2. Характеристики процесса диффузии при Т = 3000 К в различных направлениях: (а) среднеквадратичное отклонение, определяемое формулой (6) как функция времени; (б) отношение среднеквадратичного отклонения ко времени, дающее возможность оценки коэффициента самодиффузии.

Коэффициент диффузии определяется как

П = / г 2Ц)

1 п-12 2 I

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

На рис. 3 (а, б, в, г) приведено состояние атомной системы в четыре момента времени, и красным кружком выделено положение мобильного атома в

процессе эволюции. Хорошо видно, что атом соединяется с атомом, системой атомов и колеблется около нового положения некоторое время, и, если позволяет кинетическая энергия, отрывается от связи и продолжает движение. Этот процесс во времени представляет релаксационные колебания, представленные на рис. 2.

Заметим, что для графитовой конфигурации энергия связи в слое есть 7,8 эВ/атом, а между слоями энергия связи значительно меньше и находится в диапазоне 0,043-0,061 эВ/атом. Для кластерной графитовой конфигурации важную роль играют границы кластеров, которые не являются параллельными базовым плоскостям. Это обстоятельство влияет на определение коэффициента диффузии. Однако область кристаллитов нерегулярна, ее можно только оценивать. Оценка границ кристаллитов в процессе диффузии очень важна и является целью дальнейших исследований.

Рис. 3. Состояние атомной системы в четыре момента времени. Красным кружком выделено положение мобильного атома в процессе эволюции.

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

Проведено квантовое моделирование «из первых принципов» процесса диффузии атомов в аморфном углероде при структурном фазовом переходе. Показано, что термические эффекты приводят к образованию графитовой послойной структуры в молекулярной системе. В окрестности термически обусловленного фазового перехода процесс диффузии является анизотропным и обусловлен захватом атома в определенную молекулярную структуру и образованием новых ковалентных связей. Временной процесс увеличения среднеквадратичного отклонения отличается от Эйнштейновской зависимости из-за взаимодей-

а

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

Работа выполнена при поддержке гранта РФФИ №1301-12078 офи_м

This work was supported by the Russian Foundation for Basic Research, project no. 13-01-12078 ofi_m.

Список литературы:

1. Sebastian A., Pauza A., Rossel C., et al. Resistance switching at the nanometre scale in amorphous carbon // New Journal of Phy-sycs. 2011. Vol. 13. P. 013020.

2. Rueckers T., Kim K., Joselevich E., et al. Carbon nanotube based nonvolatile random access memory for molecular computing // Science. Vol. 289(5476). P. 94-97.

3. Meijer G. Who wins the nonvolatile memory race? // Science. 2008. Vol. 319(5870). P. 1625-1626..

4. Wuttig M., Yamada N. Phase-change materials for rewritable data storage // Nat. Mater. 2007. Vol. 6(12). P. 1004.

5. Standley B., Bao W., Zhang H., et al. Graphene-based atomic-scale switches // Nano Lett. 2008. Vol. 8(10). P. 3345-3349.

6. Robertson J. Diamond-like amorphous carbon // Materials Science and Engineering: R: Reports. 2002. Vol. 37(4-6). P. 129-281.

7. Silva R. Properties of amorphous carbon. The Institution of Engineering and Technology, 2003.

8. Takai K., Oga M., Sato H., et. al. Structure and electronic properties of a nongraphitic disorded carbon system and its heat-treatment effects // Phys. Rev. B. 2003. Vol. 67(21). P. 214202214212.

9. Ronning C., Griesmeier U., Gross M., et al. Conduction processes in boron- and nitrogen-doped diamond-like carbon films prepared by mass-separated ion beam deposition // Diamond and Related Materials. 1995. Vol. 4(5-6). P. 666-672.

10. Jornada F. H., Gava V., Martinotto A. L., et al. Modeling of amorphous carbon structures with arbitrary structural constraints // Journal of Physics: Condensed Matter. 2010. Vol. 22(39). P. 395402.

11. He Y., Zhang J., Guan X., et. al. Molecular Dynamics Study of the Switching Mechanism of Carbon-Based Resistive Memory // IEEE Transactions on Electron Devices. 2010. Vol. 57(12). P. 34343441.

12. Shumkin G.N., Zipoli F., Popov A.M., Curioni A. Multiscale quantum simulation of resistance switching in amorphous carbon // Procedia Computer Science. 2012. Vol. 9. P. 641-650.

13. Шумкин Г.Н., Попов А.М. Моделирование из первых принципов фазового перехода в аморфном углероде // Математическое моделирование. 2012. Т. 24, N 10. С. 65-79.

14. Попов А.М., Никишин Н.Г., Шумкин Г.Н., Многомасштабное квантовое моделирование процесса структурного фазового перехода и теплового пробоя в наноточке аморфного углерода // Computational nanotechnology. 2014. N 1. С. 17-25.

15. Marx D., Hutter J. Ab initio molecular dynamics: Theory and implementation // Modern Methods and Algorithms of Quantum

Chemistry. John von Neumann Institute for Computing, Forschungszentrum Jülich, 2000. P. 329-477.

16. Andreoni W., Curioni A. New Advances in Chemistry and Materials Science with CPMD and Parallel Computing // Parallel Computing. 2000. Vol. 26(7-8). P. 819-842.

17. Kohn W. Density Functional and Density Matrix Method Scaling Linearly with the Number of Atoms // Phys.Rev.Lett.1996. Vol. 76(17). P. 3168-3171