CC BY
54
1.3. МНОГОМАСШТАБНОЕ КВАНТОВОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА СТРУКТУРНОГО ФАЗОВОГО ПЕРЕХОДА И ТЕПЛОВОГО ПРОБОЯ В НАНОТОЧКЕ АМОРФНОГО УГЛЕРОДА
Попов Александр Михайлович, профессор, д-р физ.-мат. наук, факультет Вычислительной математики и кибернетики. Московский Государственный Университет им. М.В. Ломоносова, Москва. E-mail: popov@cs.msu.su
Никишин Николай Глебович, аспирант, факультет Вычислительной математики и кибернетики, Московский Государственный Университет им. М. В. Ломоносова, Москва. E-mail: mail@nikishinng.ru
Шумкин Георгий Николаевич, разработчик ПО, канд. физ.-мат. наук, IBM Systems and Technology Group, Научно-технический центр IBM в России, Москва. E-mail: georgiy-sh@yandex.ru
Аннотация: В научно-исследовательской лаборатории IBM в Цюрихе была проведена серия экспериментов, в которых была показана возможность создания памяти на фазовых переходах на основе аморфного углерода.
В работе предложена многомасштабная модель памяти на фазовом переходе. Процесс моделируется самосогласованно на трех уровнях. На первом уровне проводятся вычисления квантовой молекулярной динамики из первых принципов c учетом пространственного распределения температуры. Временная эволюция электронной структуры моделируется на втором уровне с использованием редуцированной молекулярной динамики Эренфе-ста в окрестности фазового перехода второго рода. На третьем уровне вычисляется новое пространственное распределение температуры с помощью уравнения теплопроводности в сплошной среде. Вычисления проводились на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ.
В работе показано, что образование графитовой послойной структуры под действием температуры приводит к локализованной по пространству зависимости проводимости. А тепловая неустойчивость, поддерживающая структуру, вызывается пространственно локализованным джоулевым источником, пропорциональным температуре. Такое поведение может объяснить возникновение s-образной формы вольт-амперной характеристики, наблюдаемой в эксперименте.
Ключевые слова: многомасштабные квантово-механические коды молекулярной динамики, фазовый переход в аморфном углероде, память на фазовых переходах, нанотехнологии, суперкомпьютер IBM Blue Gene/P
1.3. MULTISCALE QUANTUM SIMULATION OF STRUCTURE PHASE CHANGE AND THERMAL DISRUPTION IN NANODOT OF 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. E-mail: popov@cs.msu.su
Nikishin Nikolai G., postgraduate, Computational Mathematics and Cybernetics Department, Lomonosov Moscow State University, Moscow. E-mail: mail@nikishinng.ru
Shumkin Georgiy N., PhD in Physics and Mathematics, Software Developer, IBM Systems and Technology Group, IBM Science and Technology Center in Russia, Moscow. E-mail: georgiy-sh@yandex.ru
Abstract: There was a number of experiments, which showed a possibility of phase-change memory building based on amorphous carbon, carried out in IBM Zurich research laboratory.
We suggest a multiscale model of phase-change memory. A phase transition is self-consistently simulated on three different time-space levels. On the first level, we use ab initio quantum molecular dynamics calculations with taking into account temperature distribution. On the second level, time dependent evolution of the electronic density is simulated on basis of reduced Ehrenfest molecular dynamics near the line of the phase transition of the second kind. On the third level, we use a heat conduction equation in continuous media to calculate new temperature distribution. For calculations, we used the IBM BlueGene/P supercomputer installed at the Faculty of Computational Mathematics and Cybernetics of the Moscow State University.
In this paper we point, that an appearance of a graphitic layer structure from an amorphous state under the influence of temperature effects leads to a localization in space of the electric conductivity. In addition, the temperature profiles,
that maintain the structure, become unstable due to the impact of a space-localized heat source. Such a behavior could explain the appearance of s-shaped volt-ampere characteristic in a conducting nanodot during the experiment
Index terms: multiscale quantum molecular dynamics codes, phase transition in amorphous carbon, phase change memory, nanotechnology, supercomputer IBM BlueGene/P
1. ВВЕДЕНИЕ.
В области наноэлектроники при создании новых технологий используется эффект переключения молекулярной структуры [1], [2], [3], [4]. Изменение атомной структуры материала, приводящее к переключению проводимости под дествием внешних полей и температуры, лежит в основе создания памяти на фазовых переходах [5], [6], [7].
Особый интерес в качестве такого материала представляет аморфный углерод (a-C), так как допускает образование клеточной структуры проводящих наноточек. В научно-исследовательской лаборатории IBM в Цюрихе была проведена серия экспериментов [8], в которых была показана возможность создания памяти на фазовых переходах на основе a-C.
Численному моделированию структур аморфного углерода посвящен ряд работ. В работах [9], [10] проводилось моделирование на основе классической молекулярной динамики с эмпирическими потенциалами для исследования различных структур a-C. Основным недостатоком такого моделирования является многозначность появляющихся структур, что обусловлено неопределенностью в задании эмпирического потенциала взаимодействия между частицами. Квантовомеханическое моделирование применялось для изучения структур тетраэдрального, алмазоподобного аморфного углерода. Совместная модель связанных уравнений классической молекулярной динамики и транспортных уравнений сплошной среды предложена в [11] для описания памяти на фазовых переходах в аморфном углероде. Однако требуется дальнейшее исследование вопросов, связанных с согласованием поведения молекулярной системы на основе квантовой молекулярной динамики с пространственным распределением температуры, получаемой из уравнений сплошной среды.
В работах [12] и [13] проведено компьютерное квантовомеханическое моделирование процессов возникновения анизотропной атомной структуры в aC. В них показано, что появление графитовой структуры из аморфного состояния обусловлено термическими эффектами. Однако были использованы периодические граничные условия, и вычисления проводились при однородной по пространству температуре. При молекулярном моделировании движения системы молекул температура является внешним параметром. Повышение проводимости системы, связанное с
изменением молекулярной структуры в эксперименте, должно приводить к существенному увеличению джоулева нагрева проводящей точки и дальнейшему повышению температуры. Описанный процесс может привести к термическому срыву в полупроводнике и возникновению отрицательного сопротивления наноточки. Такая модель переключательного процесса может объяснить s-образное поведение вольт-амперной характеристики, наблюдаемой в экспериментах по созданию памяти на фазовых переходах. Пространственное поведение температуры должно определятся другой моделью -моделью сплошной среды. Исследования показывают необходимость многомасштабного моделирования переключения транспортных свойств атомной конфигурации.
В данной работе предлагается многомасштабная модель, связывающая квантовую молекулярную динамику системы и модель теплопроводности сплошной среды, для самосогласованного определения неоднородной по пространству температуры.
В многомасштабном моделировании на первом уровне (уровне I) для моделирования атомной структуры используется квантово-механическая молекулярная динамика Кара-Парринелло, представленная в численном коде CPMD [14], [15], основанная на теории функционала плотности [16]. Рассматривается молекулярная система,
представляющая аморфную структуру а-С. Изучается влияние температурных эффектов на формирование атомной анизотропной графитоподобной структуры в процессе молекулярной динамики. Распределенная по пространству температура в системе атомов поддерживается при помощи взаимодействия с термостатом Нозе-Хувера [14]. Основное внимание уделяется возникновению анизотропных параметров, характеризующих структуру молекулярной системы. На втором уровне (II) изучается влияние эволюции ионной структуры и внешнего электрического поля на динамику электронов. Для такого исследования используется модифицированная молекулярная динамика Эренфеста [4], [9]. На уровне (III) решается трехмерная граничная задача для уравнения теплопроводности с источником джоулева тепла, имеющим структуру, диктуемую моделями первых двух уровней.
Цель работы состоит в том, чтобы показать, что термический эффект образования молекулярной структуры из аморфного состояния приводит к
стратификации области и возникновению филаментации тока. ^ратификация атомной структуры под действием внешнего напряжения вызывает рост джоулева тепла и температурный срыв, заключающийся в дальнейшем увеличении температуры и проводимости графитоподобных слоев. Пространственно локализованный джоулев источник инициирует возникновение тепловой неустойчивости и таким образом поддерживает структуру и может объяснить форму вольт-амперной характеристики проводящей наноточки в эксперименте, проведенном в [8].
Расчеты в рамках предложенной модели требуют большой вычислительной мощности. Вычисления проводились на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ.
2. МОДЕЛИРОВАНИЕ ПЕРЕСТРОЙКИ АТОМНОЙ КОНФИГУРАЦИИ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ.
При построении моделей молекулярной динамики из «первых принципов» исходят из уравнения Шредингера и кулоновского взаимодействия заряженных частиц. Мы используем математическую модель молекулярной динамики Кар-Парринелло [14], [15], [16] на основе приближения Борна-Оппенгеймера. В приближении Борна-Оппенгеймера считается, что движение массивных ядер может быть описано в классическом приближении. Силы, действующие на ядра со стороны электронов, рассчитываются на основе решения уравнения Шредингера для электронов в основном состоянии. Мы будем использовать квантово-механическую формулировку, основанную на уравнениях функционала плотности («DFT» [16]).
Уравнения движения ядер в усредненном поле электронов записываются в виде:
m, RI (0 = -£V
Z,Z,
Z,
R i - r |
dr '
(1)
(" 2 A + U eff Г ) К (r) = SnVn (r) ' Плотность электронов записывается как:
Ne
p(r) = XiKn (r)|2,
(2)
(3)
и ,ff(p, r) = Uh (p,r) + U xc(p,Vp,r)+ 2
—
r)- 2 -■-+ ,
71 i Ri - r i i R, - R, i
, (4)
где ин - потенциал Хартри, а ихс — обменно-
корреляционный потенциал в форме PBE (Perdew-Burke-Ernzerhof) [14].
Для поиска электронных орбиталей уп возникает
необходимость решения нелинейной задачи на собственные значения. Эффективный потенциал иеП(р,г) зависит от электронной плотности р, которая, в свою очередь, зависит электронных орбиталей уп. При известном начальном приближении для потенциала по решению задачи на собственные значения еп и уп из формулы (3) может быть найдена пространственная зависимость электронной плотности системы. А новая пространственная зависимость эффективного потенциала и^р,г) может быть определена по рассчитанной плотности.
Электронные орбитали уп раскладываются по плоским волнам в приближении нулевой зоны Бриллю-
эна:
Wn (r) = ^=Z cn,k exp(/k •r ),
(5)
Для поиска новых орбиталей требуется решать уравнение Шредингера при каждом новом положении ядер:
где величина Ме равна количеству валентных электронов в системе. Эффективный потенциал представляется в форме:
где Q - объем ячейки, k — вектор обратной решетки, и cn,k — Фурье-коэффициенты орбитали n.
3. ЧИСЛЕННАЯ ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ КВАНТОВОЙ МОЛЕКУЛЯРНОЙ ДИНАМИКИ НА СУПЕРКОМПЬЮТЕРЕ IBM BLUE GENE/P.
Для расчетов фазового перехода при помощи кода CPMD [15] использовался суперкомпьютер IBM Blue Gene/P, установленный на факультете ВМК МГУ. Самой затратной по времени вычислений частью является параллельный расчет электронной структуры системы. Используемый алгоритм обладает вычислительной сложностью, масштабируемой как o(N ).
Процесс решения задачи включает в себя нескольких этапов. Сначала по коэффициентам cnk при помощи трехмерного быстрого преобразования Фурье вычисляются элементы матрицы гамильтониана H. Затем решается нелинейная задача на собственные значения, и проводится пересчет плотности и эффективного потенциала в реальном пространстве с использованием трехмерного быстрого преобразования Фурье. Для поиска собственных значений на каждой итерации в коде CPMD реализован параллельный алгоритм Ланцоша с использованием технологии MPI. В данной реализации используется быстрая раздельная сеть коллективных межпроцессорных сообщений суперкомпьютера IBM Blue Gene/P.
Алгоритм решения задачи осуществляет распараллеливание на нескольких уровнях. На распределенной памяти используется крупноблочное распараллеливание. Процедуры MPI распределяют коэффициен-
n=l
ты Фурье волновых функций всех электронных состояний по узлам суперкомпьютера. Внутри крупных блоков происходит распараллеливание длинных циклов на общей памяти узла при помощи технологии OpenMP. При этом данные распределяются по межпроцессорным соединениям трехмерного тора IBM Blue Gene/P оптимальным образом.
Коэффициенты Фурье волновых функций распределяются блоками по одному из двух индексов, по k или n. По индексу вектора обратной решетки k проводится распределение базисного набора плоских волн. По орбитальному индексу n производится вычисление Фурье-преобразования волновой функции на одном узле. В реальном пространстве производится прямая декомпозиция по слоям.
Размер суперъячейки, в которой располагались атомы углерода, составлял 14.268 х 10.701х 10.701 Ä.
Начальная конфигурация включала 192 атома: 44
2
атома sp с тремя связями, представлявшие графитовые области; 78 атомов sp3 с четырьмя связями, соответствовавшие алмазным областям; 70 атомов sp с двумя связями, представлявшие линейные области. Соответствующее процентное содержание числа атомов в начальной конфигурации составляло 23%, 41% и 36%. Плотность системы была равна n = 2.34 г/см3. При этом использовалось разложение по 384 электронным орбиталям, 54554 плоских волн для каждой волновой функции и 436234 плоских волн для плотности. Временной шаг был равен 0.125 фс. Пространственная сетка состояла из 144 х 108 х 108 точек. Длина производственного цикла составляла 150 пс/неделя.
Ряд расчетов молекулярной динамики при разных температурах был проведен в работе [12]. Во всех расчетах при поднятии температуры азмазные области замещались алмазных графитовыми. Сравним установившиеся атомные структуры в процессе молекулярной динамики.
При температурах порядка T = 3000 K формируется графитоподобная послойная структура, обладающая меньшим сопротивлением, чем начальная конфигурация.
Эти изменения, обусловленные температурным эффектом, рассматриваются как фазовый переход в наноструктуре материала, приводящий к пороговой смене сопротивления. В ходе нагрева обнаруживается порог по температуре, выше которого возникают высокопроводящие кластеры графита.
Две конфигурации, соответствующие расчетам молекулярной динамики, приведены для температуры T = 3000 K на рис. 1 (а) в начале нагрева, когда конфигурация является аморфной, и на рис. 1 (б) в конце процесса, когда образуются графитоподобные структуры. Возникновение высокопроводящих атомных конфигураций
происходит в результате замещении алмазных областей графитовыми.
Для перехода в высокопроводящее состояние, как следует из расчетов, необходим не только рост числа графитоподобных атомов, но и то, чтобы атомы соединялись в графитовые слои с кольцами атомов углерода на них. В переключении сопротивления важной характеристикой является именно рост графитовых кластеров.
Структура, полученная в конце расчета молекулярной динамики, в среднем отражает параметры структуры графита: энергия связи между слоями равна 16.8 Дж/моль, расстояние между слоями равно 0.3 нм, энергия связи в слое 167.6 Дж/моль (1118 °С), расстояние между атомами равно 0.1 нм.
Введем параметры структуры, характеризующие образование порядка из аморфной молекулярной структуры. Для этого рассмотрим пространственную плотность ионов в процессе молекулярной динамики. Пространственная плотность ионов вычислялась с использованием техники гистограмм в процессе эволюции системы. Направление вдоль оси х выберем как перпендикулярное графитовым слоям, а направления у и z как параллельные графитовым плоскостям. Гистограммы плотности Nа получаются
как отношение количества атомов, находящихся внутри пространственных слоёв, образованных параллельными плоскости Oyz плоскостями, к общему числу атомов в системе. Гистограмма вдоль оси х аппроксимируется гауссианами с центрами в точках, соответствующих положению графитовых плоскостей. На рис. 2 показана плотность ионов, характеризующая эволюцию атомной конфигурации при Т = 3000 К, в начале процесса нагрева, а на на рис. 3 - в конце.
Рис. 1. Эволюция структуры молекулярной системы при температуре Т = 3000 К (а) аморфная конфигурация в начале процесса нагрева, (б) графитоподобная конфигурация в конце процесса.
Рис. 2. Плотность ионов для аморфной конфигурация в начале процесса нагрева, полученная методом гистограмм. Распределение плотности аппроксимировано суммой гауссиан.
На основе статистики из вычислений СРМй, используя положения атомов при определенной температуре, рассчитывается структурный фактор 8(/). В течение эволюции структура локализуется в окресности графитовых слоев, что приводит к филаментации тока через образец. В качестве структурного параметра можно выбрать полуширину гауссиан, аппроксимирующих гистограммы плотности. Временная эволюция такого структурного фактора 8 (/) показана на рис. 4.
Рис. 3. Плотность ионов для графитоподобной конфигурации в конце процесса, полученная методом гистограмм. Распределение плотности аппроксимировано суммой гауссиан.
Рис. 4. Временная зависимость структурного фактора и её аппроксимация (показана линией)
Эволюция плотности ионов в процессе молекулярной динамики позволяет вычислить профиль свободной энергии в процессе переключения структуры.
4. РЕДУЦИРОВАННАЯ МОДЕЛЬ ЭВОЛЮЦИИ ЭЛЕКТРОННОЙ ПОДСИСТЕМЫ НА ОСНОВЕ МОЛЕКУЛЯРНОЙ ДИНАМИКИ ЭРЕНФЕСТА.
На уровне (II) учитывается временная эволюция электронной структуры под действием приложенного напряжения. В модели этого уровня рассчитывется не полная временная эволюция, а только эволюция в окрестности структурного фазового перехода второго рода. Необходимые структурные параметры и электронные орбитали берутся из модели первого уровня. Введение внешнего электрического
потенциала приложенного напряжения позволяет сделать оценки электрического тока и его зависимости от величины приложенного напряжения и критериев переключения. Это может быть сделано в рамках модифицированной молекулярной динамики Эренфеста [4] , [12], [13].
Рассмотрим конструкцию редуцированной модели молекулярной динамики Эренфеста. Электронные волновые функции раскладываются по полному набору орбиталей:
Т (Г, R, t ) = Yfm (t )Vm (Г, R)
(6)
где у/т (г, R) удовлетворяет стационарному уравнению:
И¥т (Г, R) = ЕпУт (Г, R) . (7)
Молекулярная динамика Эренфеста имеет дело с временными зависимостями коэффициентов Ст) .
Мы вводим структурный параметр s(t), который соответствует положению ионов при заданной температуре. Рассмотрим теперь уравнение для эволюции ионов и электронов в непосредственной окрестности структурного фазового перехода. Приближенные уравнения движения можно записать в следующей форме для структурного параметра s(t) :
d s
dF
dE
, (8)
^ = "ЕК -^рп СК ШЕп - Ек )d пк (8)
и для эволюции электронов:
¡М = ^ (ОЕ (8) - ^Ск(^| • Лпк(0' (9)
где г есть положения электронов, а положения ядер К1 участвуют в уравнении через структурную переменную ). Уравнения включают следующие интегралы перекрытия электронных орбиталей:
д
dL = JVm dXc^kdr
dlk = Jvm jyVkdr' dík = JVm ~^Vkdr
vex=V i vext i Vk) •
(10)
(11)
Связь электронных и ионных уравнений осуществляется за счет интегралов перекрытия орбиталей dтк, а каждое состояние характеризуется
своей энергией Кона-Шэма, которая соответствует определенной парциальной частоте колебаний. Связь осцилляторов определяется интегралом перекрытия Утк, связанным с внешним потенциалом. Внешний потенциал зададим в следующей форме: ^ = V0 • (2 - ^о)ехр(-(г - Г0)2/Ст02). (12)
Перемещения ионов через структурный параметр учитывается вторым членом справа в электронных уравнениях. При отсутствии перемещения электронные уравнения представляют собой линейные уравнения
связанных осцилляторов. Полная система при этом оказывается нелинейной за счет члена, связанного с изменением структурного параметра.
Рассмотрим двухмодовую систему, включающую основное состояние, характеризующееся орбиталью HOMO (Highest Occupied Molecular Orbital, высшая занятая молекулярная орбиталь) и энергией E1, и возбужденное состояние с орбиталью LUMO (Lowest Unoccupied Molecular Orbital, низшая незанятая молекулярная орбиталь) и энергией E2 . Такая система уравнений принимает наиболее наглядную форму. Введем следующие переменные.
_ Л „2
Плотности электронных состояний р = Cj + CR и
22
Р2 = C2 + CR , причем переключение системы из одного состояния в другое будем характеризовать разностью u = р -р2;
Вспомогательную переменную ß = CRCR + C/C2;
Вспомогательную переменную у = CRCR - CRCR;
Ширину запрещенной зоны G = E2 -Ej;
Интеграл перекрытия, связанный с внешним приложенным потенциалом, обозначим через
Также введем вспомогательную переменную В этом случае система уравнений примет вид:
du dt
4V ■/ + 4 f-ß,
dß r , — = gy-f ■u, dt
^ = -G-ß-V■ u ■ dt
(13)
У получившейся системы есть интеграл движения:
2u2 + 2ß2 + 2y2 = const.
(14)
Переменная у связана с полным током, протекающим через сечение полупроводника: I = -у d . (15)
На рис. 5 изображен процесс переключения и возникновения тока при быстром изменении параметра структуры со скоростью, имеющей форму гауссиана с полушириной равной 0.5. Кривая 1 соответствует переменной и, кривая 2 - переменной в, кривая 3 - переменной у, кривая 4 - параметру f (скорость изменения структуры) и 5 - интегралу движения. В этом случае происходит полная перекачка энергии из основного состояния в возбужденное.
m
V = V = V ; v ' 1 2 V 2 1 '
11(1)
-
' 1
1 \ ^
- \ ^^2 ■
О 5 10 15 20
I
Рис. 5. Процесс переключения для двух состояний при скорости изменения структуры, имеющей форму гауссиана с полушириной равной 0.5. Кривая 1 соответствует переменной и, кривая 2 - переменной в, кривая 3 - переменной у, кривая 4 - параметру f (скорость изменения структуры) и 5 - интегралу движения.
На рис. 6 представлено переключение при продолжительном росте структуры, который апрроксимирован квадратичной зависимостью по времени. Кривая 1 соответствует переменной и, кривая 2 - переменной в,кривая 3 - переменной у, кривая 4 - параметру f (скорость изменения структуры) и 5 - интегралу движения. Ток выходит на стационарное значение при наличии быстрых осцилляций в электронной плотности.
и(1)
О 5 10 15 20
I
Рис. 6. Процесс переключения при продолжительном росте структуры, аппроксимированной квадратичной зависимостью по времени. Кривая 1 соответствует переменной и, кривая 2 - переменной в, кривая 3 -переменной у, кривая 4 - параметру f (скорость изменения структуры) и 5 - интегралу движения. В процессе переключения коэффициенты системы электронных уравнений меняются в соответствии с анизотропным характером переключения. Коэффициенты уменьшаются поперек базовых
поверхностей и становятся на порядок меньше продольных коэффициентов.
5. МАКРОСКОПИЧЕСКАЯ МОДЕЛЬ ТЕПЛОПРОВОДНОСТИ В СПЛОШНОЙ СРЕДЕ.
Температурный эффект обуславливает структурный фазовый переход в молекулярной структуре. Но для модели молекулярной динамики температура является заданной термодинамической постоянной, как по времени, так и по пространству. Это в макроскопическом смысле противоречит эксперименту, где фазовый переход происходит в проводящей точке диаметром порядка 20 нм. Пространственная зависимость температуры и учет границ проводящей точки могут существенно влиять на результаты моделирования. Температура для расчетов молекулярной динамики, по нашим предположениям, должна диктоваться
макроскопическим уровнем модели с учетом граничных условий. Для моделирования на макроскопическом уровне в работе используется уравнение теплопроводности с источником джоулева тепла, генерируемого плотностью тока ],
напряженностью поля Е и пространственным распределением проводимости. Необходимо согласование распределения температуры в модели теплопроводности и передача ее в модель молекулярной динамики. Такое моделирование проводится для окрестности только одного графитового слоя. Полученная температура затем передается в молекулярную динамику для необходимого числа слоев. В то же время структурные параметры, определяющие
электропроводность, должны быть переданы в макроскопическую модель теплопроводности. Эта процедура повторяется через определенные промежутки времени. Уравнение теплопроводности принимает следующий вид:
рС д-Г = 1V1T) + V,, (K||V||r) + q0 (] • Е), (166)
где р - плотность, С - теплоемкость, К± -коэффициент теплопроводности в направлении М -перпендикулярном графитовым слоям, К, -коэффициент теплопроводности вдоль графитовых слоев ^,у), q0 - безразмерный параметр. Закон Ома
записывается в виде произведения тензора проводимости Г (продольная и поперечная составляющие) на напряженность поля Е.
] = г(Е,Т, з(? ))• Е. (177)
Тензор проводимости г(Е,Т, учитывает
пространственную зависимость, зависимость от температуры и структурного фактора ), полученных при моделировании на микроскопическом уровне.
Из результатов расчетов следует наличие эффекта локализации тепла в определенной конечной области. Зависимость электропроводности от температуры и пространственных переменных определяет границы по параметрам мощности и структурного параметра, при которых наступает температурная неустойчивость, и начинается неограниченный рост температуры в центре ячейки. Этот эффект мы назовем температурным пробоем или температурным срывом. Механизм неустойчивости состоит в том, что вкладываемая источником мощность становится больше, чем диссипация тепла за счет теплопроводности. При линейной зависимости проводимости от температуры порог срыва определяется как граница температурной неустойчивости.
Порог по 8 равен 8ш =0.1. Это значение найдено
в молекулярных расчетах в конце нагрева. Профили температуры в поперечном к гафитовым слоям х-направлении и продольном z-направлении при 8 = 0.1 показаны на рис. 7 и рис. 8.
Рис. 7. Профили температуры в поперечном направлении (x-направлении) при S = 0.1.
Рис. 8. Профили температуры в продольном направлении (z-направлении) при S = 0.1 .
Рис. 9. Молекулярная структура при распределенной по пространству температурой, полученной в модели теплопроводности. Алмазоподобные слои соответствуют границам проводящих графитоподобных слоев.
Пространственные зависимости, полученные в модели теплопроводности, учитывались в квантовой молекулярной модели. На рис. 9 показана образующаяся молекулярная структура, при температуре в центре порядка T = 4500 K и на границе порядка T = 600 K. Для проведения расчетов с распределенной температурой использовалось взаимодействие с термостатом Нозе-Хувера [14]. Главный эффект состоит в том, что в центральной области, области с увеличивающейся температурой, конфигурация атомом в перестраивается в графитоподобную, что приводит к повышению проводимости в центре. В свою очередь, на границе проводящей точки преобладает алмазоподобная структура.
6. ЗАКЛЮЧЕНИЕ.
В работе предложена многомасштабная модель памяти на фазовом переходе и проведено моделирование из первых принципов перехода молекулярной структуры аморфного углерода в графитоподобную конфигурацию под действием температуры, направленное на объяснение механизма смены сопротивления в экспериментах IBM [8].
В результате эволюции электронной подсистемы во внешнем электрическом поле в системе возникнает ток, что соответствует возрастанию продольной проводимости графитоподобной структуры. Перестройку атмоной конфигурации можно рассматривать как фазовый переход второго рода в наноструктуре материала, приводящий к пороговой смене электрической проводимости.
Для сравнения с экспериментом существенно необходим учет пространственного распределения температуры. В работе показано, что образование графитовой послойной структуры под действием температуры приводит к локализованной по пространству зависимости проводимости. А тепловая неустойчивость, поддерживающая структуру,
вызывается пространственно локализованным джоулевым источником, пропорциональным температуре. Такое поведение может объяснить возникновение s-образной формы вольт-амперной характеристи, наблюдаемой в эксперименте.
Работа выполнена при поддержке гранта РФФИ №1301-12078 офи_м
This work was supported by the Russian Foundation for Basic Research, project no. 13-01-12078 ofi_m
Список литературы:
[1] Rieth M., Schommers W. Handbook of Theoretical and Computational Nanotechnology. Karlsruhe, Germany: Forschungszentrum Karlsruhe, 2006.
[2] Попов А.М. Вычислительные нанотехнологии. М. : КНОРУС, 2014.-312c.
[3] Liljeroth P., Molen S. Charge transport through molecular switches // Journal of Physics: Condensed Matter. 2010. 22, N 13, 133001-133030
[4] Shumkin G.N., Popov A.M., Curioni A., Laino T. A multiscale modeling of naphthalocyanine-based molecular switch // Procedia Computer Science, 2010, v. 1, no. 1, p. 185-192.
[5] Wuttig M., Yamada N. Phase-change materials for rewritable data storage // Nat. Mater. 2007. 6. 824-832
[6] Standley B., Bao W., Zhang H., et al. Graphene-based atomic-scale switches // Nano Lett. 2008. 8, N 10. 3345-3349
[7] Meijer G. Who wins the nonvolatile memory race? // Science. 2008. 319. 1625-1626
[8] Sebastian A., Pauza A., Rossel C., Shelby R.M., Rodriguez A.F., Pozidis H., Eleftheriou E. Resistance switching at the nanometre scale in amorphous carbon // New Journal of Physycs. 2011. 13. 013020
[9] 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. 67. 214202214212
[10] Jornada F. H., Gava V., Martinotto A. L., Cassol L. A., Perot-toni C. A. Modeling of amorphous carbon structures with arbitrary structural constraints // Journal of Physics: Condensed Matter. 2010. 22, N 39. 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. 57, N 12. 3434-3441
[12] Shumkin G.N., Zipoli F., Popov A.M., Curioni A. Multiscale quantum simulation of resistance switching in amorphous carbon // Procedia Computer Science. 2012. 9. 641-650.
[13] Шумкин Г.Н., Попов А.М. Моделирование из первых принципов фазового перехода в аморфном углероде // Математическое моделирование. 2012. Т.24, №10. 65-79.
[14] Marx D., Hutter J. Ab initio molecular dynamics: Theory and implementation // Modern Methods and Algorithms of Quantum Chemistry. 2000. 1. 329-477
[15] Andreoni W., Curioni A. New Advances in Chemistry and Materials Science with CPMD and Parallel Computing // Parallel Computing. 2000. 26. 819-842
[16] Kohn W., Density Functional and Density Matrix Method Sacaling Linearly with the Number of Atoms // Phys.Rev.Lett.. 1996. 76. 3168-3171
