CC BY
16
УДК 517.958:518.12(075.8)
А.М. Попов, Г.Н. Шумкин2
МОДЕЛИРОВАНИЕ ПЕРЕКЛЮЧЕНИЯ СОПРОТИВЛЕНИЯ В АМОРФНОМ УГЛЕРОДЕ*
На основе квантовой молекулярной динамики проведено математическое моделирование механизма переключения сопротивления в аморфном углероде (а-С). Электрическая проводимость в тонких пленках аморфного углерода связывается с фазовым переходом, обусловленным электрическим полем и сопутствующими джоулевым нагревом и охлаждением. Моделирование показывает, что механизм перехода связан с кластеризацией существующих проводящих областей графита (sp2) внутри непроводящей алмазной (sp3) матрицы. Из моделирования найдены условия для фазового перехода в а-С из одного состояния в другое. Были также рассчитаны возбужденные молекулярные орбитали и ширина электронной щели при различных температурах.
Ключевые слова: аморфный углерод, материал фазового перехода, квантовая молекулярная динамика, переключение сопротивления, переключение проводимости.
1. Введение. В последнее время много внимания уделяется использованию различных форм углерода, таких, как углеродные нанотрубки и графен, в создании энергонезависимой памяти [1]. По сравнению с этими формами углерода в особенности интересным является аморфный углерод (а-С) как наиболее простой, недорогой и хорошо известный материал для создания энергонезависимой памяти, основанной на фазовых переходах [2]. В последних экспериментах, проведенных в научно-исследовательской лаборатории IBM в Цюрихе [3], найдено подобие между поведением тонких пленок а-С и поведением памяти на фазовых переходах из традиционных материалов, таких, как халькоге-нидные стекла (GST). К халькогенидным стеклам, например, относятся СегЭЬгТеб, GeTe и другие широко исследуемые сплавы из таких же элементов. Цель настоящей работы — моделирование из первых принципов экспериментов по переключению сопротивления в аморфном углероде, проведенных в [3].
Уникальные свойства аморфного углерода связываются с особенностями его мезоскопического строения [4]. Пленки аморфного углерода состоят из двух встроенных друг в друга фаз — алма-зоподобной, характеризующейся «р3-гибридизацией атомов углерода и тетраэдрической структурой, которая обладает малой электропроводностью, и графитоподобной фазой, состоящей из фрагментов графитовых плоскостей и изогнутых фуллереноподобных фрагментов. Графитоподобная фаза характеризуется «р2-гибридизацией углерода и обладает сравнительно высокой электропроводностью в зависимости от размеров кластеров. Фазы перемешаны фрагментами с характерным размером от единиц до сотен ангстрем.
Найдено, что изменение отношения sp3/sp2 сильно влияет на проводимость углерода. Но даже при одном и том же отношении sp3/sp2 кластеризация углерода sp2 может приводить к значительному увеличению проводимости. Определению размеров кластеров в такой структуре посвящено немало работ, например [4]. Подчеркивается, что свойства графитоподобных кластеров очень сильно зависят от размеров фрагментов.
Численному моделированию структур аморфного углерода посвящено много работ. Моделирование на основе классической молекулярной динамики с эмпирическими потенциалами (например, Тер-соффа) использовалось для описания различных структур а-С [5]. Моделирование из первых принципов применялось для изучения структур тетраэдрического аморфного углерода [6]. Совместная модель связанных уравнений классической молекулярной динамики и транспортных уравнений сплошной среды предложена в [7] для описания памяти на фазовых переходах в аморфном углероде. Тем не менее такой механизм переключения, как кластеризация, не понят до конца.
1 Факультет ВМК МГУ, проф., д.ф.-м.н., e-mail: popovQcs.msu.su
2 Факультет ВМК МГУ, асп., e-mail: georgiy-shQyandex.ru
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проект № 11-01-00216-а.
Исследования, представленные в настоящей работе, направлены на суперкомпьютерное математическое моделирование процесса фазового перехода в аморфном углероде. Моделирование основано на вычислениях из первых принципов на основе квантово-механической молекулярной динамики Кар-Парринелло, реализованной в численном коде CPMD [8, 9].
Вычисления проводились на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ.
2. Постановка задачи. При построении моделей "молекулярной динамики из первых принципов" исходят из уравнения Шредингера и кулоновского взаимодействия заряженных частиц. Мы используем математическую модель молекулярной динамики Кар-Парринелло [10] на основе приближения Борна-Оппенгеймера. В приближении Борна-Оппенгеймера считается, что движение массивных ядер может быть описано в классическом приближении. Силы, действующие на ядра со стороны электронов, рассчитываются на основе решения уравнения Шредингера для электронов в основном состоянии. Мы будем использовать квантово-механическую формулировку, основанную на уравнениях функционала плотности (БИ1) [11]. Уравнения модели Кар-Парринелло записываются в виде
д
М/К/ф = [{-фп}, {И/}],
ß4>n(t) =
6
5-фп
-EKS [Ш, {R-/}] + У^ КтФ
т
(1)
где М/ и К/ — масса и координата ядра I, ц — фиктивная масса электрона, фп — орбитали Кона-Шэма и Е'кв — энергия Кона-Шэма. Последний член в уравнении (1) есть голономные связи, которые учитываются требованием ортогональности орбиталей. Это приводит к необходимости вычислять си-
ТШ^Екв)-, орбитали {^-¡щргЕкв ) и множители Лагранжа Апт на каждом
лы, действующие на ионы шаге моделирования.
Энергия Кона-Шэма определяется выражением
Ne
где L^xt есть ионный потенциал и
dr+
Uext(r)p(r) dr + 1 ff Р(ГГ)РУ drdr' + Exc[p] + ВДН,}],
Ne
(2)
есть плотность заряда электронов, Ые — число электронов, /п — число заселенности, Ехс — обменно-корреляционная энергия и /•,'/ — электростатическое ион-ионное взаимодействие.
3. Численная параллельная реализация квантовой молекулярной динамики на суперЭВМ IBM Blue Gene/P. Расчет фазового перехода реализован на суперЭВМ IBM Blue Gene/P, установленной на факультете ВМК МГУ. Наиболее затратной в смысле времени вычислений является часть задачи, связанная с параллельным расчетом электронной структуры системы. В рамках теории функционала плотности (DFT) вычисления сводятся к следующим этапам. Уравнения Кона-Шэма приводят к необходимости решения нелинейной задачи на собственные значения для волновых функций:
Нфп( г) = еп-фп( г),
где гамильтониан равен
Я=- QA + £/eff(p,r)) •
Нелинейность задачи связана с тем, что эффективный потенциал Ueg(p, г) зависит от электронной плотности р, которая в свою очередь зависит от решения задачи — волновых функций. Поэтому необходим итерационный процесс. Допустим, что на начальной итерации известно некоторое приближение
для потенциала. Тогда, используя решение задачи на собственные значения вп и собственные функции фп, находим пространственную зависимость электронной плотности системы по формуле (2). По рассчитанной плотности определяется новая пространственная зависимость эффективного потенциала
Вычислительная сложность используемого алгоритма масштабируется как 0(N3), где Ne — число электронов в системе. Для нахождения собственных значений на каждой итерации в коде CPMD реализован параллельный алгоритм Ланшоца с использованием технологии MPI [12]. Его особенностью является использование быстрой раздельной сети для коллективных сообщений суперкомпьютера IBM Blue Gene/P.
Для решения задачи в коде CPMD используется быстрое преобразование Фурье для преобразования величин между двумя пространствами. Это связано с тем, что при представлении решения в виде ряда Фурье оператор кинетической энергии, представляющий сумму вторых производных, приводит просто к умножению на квадрат волнового числа. При этом оператор потенциальной энергии в реальном пространстве представляет просто умножение на потенциал, а в фурье-пространстве приводит к необходимости вычислять свертки. Чтобы избежать этого, вычисление эффективного потенциала и его умножение на волновую функцию проводятся в реальном пространстве.
Таким образом, рассматривая представление нулевой зоны Бриллюэна, орбитали Кона-Шэма записываются в виде
где Q — объем ячейки, к — вектор обратной решетки, сп^ фурье-коэффициенты орбитали п. Используется следующая схема решения. После интегрирования по коэффициентам сП)к с использованием 3dFFT (быстрое трехмерное преобразование Фурье) вычисляются элементы матрицы гамильтониана Н. Затем решается нелинейная задача на собственные значения. Далее происходит пересчет плотности и эффективного потенциала в реальном пространстве с использованием 3dFFT.
Несколько уровней распараллеливания определяют используемые параллельные стратегии. Крупноблочное распараллеливание используется на распределенной памяти. С помощью процедур MPI производится распределение коэффициентов Фурье волновых функций для всех электронных состояний на все процессоры. Распределение данных производится так, что минимизируется число обменов с поддержанием загрузки в обоих пространствах. Для распараллеливания длинных циклов используется ОрепМР-распараллеливание на общей памяти узла. Также данные могут быть оптимально распределены по трехмерному тору IBM Blue Gene/P.
Блочное распределение коэффициентов Фурье волновых функций проводится либо по орбитальному индексу п, либо по индексу вектора обратной решетки к. По п производится вычисление фурье-преобразования волновой функции на одном узле. Здесь имеется предел по размеру ячейки моделирования и размеру набора базисных функций из-за памяти, приходящейся на узел. По к проводится распределение базисного набора плоских волн. При этом необходимо выполнять условия поддержания баланса загрузки. В реальном пространстве производится прямая декомпозиция по слоям. Здесь число процессоров не должно превышать наибольший размер сетки для того, чтобы обеспечить баланс загрузки. Если процессоров больше, то слои дублируются по процессорам.
На IBM Blue Gene/P в МГУ при моделировании фазового перехода мы получили масштабирование системы из 192 атомов до 512 процессоров с 70%-й эффективностью и длительностью производственного цикла порядка 150 ps в неделю. При этом рассматриваются 384 орбитали, используются 54 554 плоских волны для каждой волновой функции и 436 234 плоских волны для плотности. Временной шаг равен 0.125 fs, пространственная сетка составляет 144 х 108 х 108 точек. Время вычислений для одного варианта с использованием 512 процессоров IBM Blue Gene/P порядка 24 часов.
4. Кластеризация существующих областей графита sp2 внутри алмазной матрицы sp3 как механизм переключения сопротивления. Мы исследовали возможность переключения сопротивления как термически обусловленный процесс. Для проведения процесса молекулярной динамики была выбрана начальная конфигурация атомов углерода, показанная на рис. 1, а. Начальная
к
суперъячейка состоит из области графита вр2 и области алмаза вр3. Периодическое продолжение суперъячейки представлено на рис. 1, б. Атомная конфигурация в суперъячейке состоит из 192 атомов: 44 атомов с тремя связями, обозначаемых как Ывр2 и представляющих графитовые области; 78 атомов с четырьмя связями, обозначаемых как Ывр3 и представляющих алмазные области; 70 атомов с двумя связями, обозначаемых как Ивр и представляющих линейные области. Соответствующее процентное содержание числа атомов в начальной конфигурации равно 23, 4, 36. Плотность системы равна 2.34 Размер суперъячейки составляет 14.268 х 10.701 х 10.701 А.
Рис. 1. Начальная конфигурация молекулярной системы, состоящая из графитовых вр2- и алмазных вр3-областей: а — атомная конфигурация в расчетной суперъячейке; б — конфигурация, представляющая периодическое продолжение суперъячейки
Была проведена серия расчетов молекулярной динамики системы при разных температурах. Во всех расчетах при поднятии температуры наблюдался эффект замещения алмазных областей графитовыми. Сравним установившиеся атомные структуры после процесса молекулярной динамики. Две конфигурации, соответствующие расчетам молекулярной динамики, приведены на рис. 2, а для температуры Т = 2000 К и на рис. 2, б для температуры Т = 3700 К. Видно, что при высокой температуре атомная конфигурация представлена в основном графитовыми областями.
Рис. 2. Геометрия молекулярной системы после процесса молекулярной динамики, проведенной при Т = 2000 К (а) и 3700 К (б)
Эффект замещения графитовыми областями хорошо виден из эволюции числа атомов с различными связями в процессе молекулярной динамики при температурах Т = 400 и 3700 К, что представлено на рис. 3, а. Следует заметить, что при фазовом переходе происходят два процесса на разных временных масштабах. Первый — это регибридизация атомов углерода и переход из 5]?3-гибридной конфигу-
рации в -гибридную конфигурацию. Этот процесс происходит очень быстро в течение нескольких сотен фемтосекунд. Второй процесс переноса и диффузии протекает значительно медленнее — десятки пикосекунд.
Рис. 3. Изменение числа атомов с различными связями в двух расчетах молекулярной динамики при Т = 400 и 3700 К (а). Сравнение энергии Кона-Шэма Екб молекулярной системы при различных температурах в единицах (еУ, К) (б)
Найдены два различных типа конечной атомной структуры. Первый тип структуры получен при температурах ниже Т = 2000 К. Он имеет сложную структуру, включающую смесь вр2- и .^-подобных атомов. Во всех случаях повышения температуры число 5]?2-подобных атомов возрастает, а число -подобных атомов убывает по сравнению с начальной конфигурацией. Однако вр2-атомы не образуют правильной графитовой структуры в виде слоев из колец. Часть атомов в начально выбранных графитовых слоях присоединяется к алмазоподобным областям. Конфигурации, полученные в этом диапазоне температур, при нагреве обладают большим электрическим сопротивлением. Второй тип конфигураций возникает при температурах больше Т = 2000 К. Постоянный рост числа атомов вр2 приводит к формированию графитоподобной послойной структуры. В таких структурах число атомов вр3 очень мало. Структуры второго типа обладают малым сопротивлением. Однако структура не является идеально графитовой, так как из-за наличия температуры графитовые слои искривлены, и необходима процедура отжига.
Из предыдущих расчетов видно, что изменения в кластерной структуре а-С существуют ниже точки плавления аморфного углерода. Это изменение в кластерной структуре можно рассматривать как фазовый переход в наноструктуре материала, приводящий к пороговой смене сопротивления. Сценарий с нагревом приводит к порогу по температуре, выше которого возникают высокопроводящие области. Температура является основной причиной неустойчивости кластерной структуры и приводит к появлению высокопроводящих кластеров графита внутри алмазной матрицы. В эксперименте джоулев нагрев приводит к возрастанию температуры до пороговых значений. Этот сценарий соответствует нижней части Э-образной вольт-амперной характеристики, получаемой в эксперименте [3]. Но после этого проводящая конфигурация не разрушается и существует в состоянии с высокой проводимостью. Эта стадия соответствует верхней части Э-образной вольт-амперной характеристики, в которой ток продолжает возрастать даже в условиях, когда начался процесс отжига.
Полученная немонотонная зависимость числа атомов вр2 в атомной конфигурации от температуры определяется также скоростью подъема температуры и скоростью охлаждения системы. Вычисления показывают, что для быстрого роста графитового кластера необходимо начальное присутствие определенной доли атомов вр2 в алмазоподобной структуре. Только алмазоподобные области не приводят к формированию графитовых областей. При больших температурах (больше 2000 К) количество атомов, составляющих графитовые области, возрастает, что приводит к падению сопротивления в конечном состоянии. Было получено также сравнение энергий Кона-Шэма атомных конфигураций, полученных при повышении температуры и при охлаждении (отжиг) из двух различных конфигураций. Зависимость энергии систем от температуры показана на рис. 3, б. Ясно виден ги-стерезисный характер зависимостей энергии. Это означает, что, прогревая начальную конфигурацию и затем охлаждая ее, мы не возвращаемся к исходной системе. Верхняя часть кривой соответствует повышению температуры атомной системы. Нагретая до максимальной температуры (Т = 3700 К)
атомная конфигурация затем охлаждалась. Важно отметить, что отжиг начинался с температуры выше пороговой, начиная с которой появляются графитовые структуры.
Полученные атомные конфигурации при охлаждении соответствуют графитовым состояниям. Сравнение атомных структур при одинаковой температуре (например, при Т = 2000 К), но полученных с помощью нагрева и охлаждения, показывает, что на верхней ветви расположены аморфные непроводящие конфигурации, а на нижней ветви — графитоподобные проводящие конфигурации.
Интересно получить противоположный сценарий перехода от жидкого состояния к твердому аморфному состоянию. Переход от состояний с высокой проводимостью к состояниям с низкой проводимостью проводился мгновенным охлаждением пленки аморфного углерода, полученной при высокой температуре. Этот сценарий предложен в работе [6] для аморфного алмаза. Проводилась термализа-ция при температуре Т = 5000 К начиная с конфигурации, полученной при Т = 3700 К. Уже в этом процессе число атомов sp2 сокращается, а число атомов sp3 возрастает. Затем жидкость мгновенно охлаждалась до температуры Т = 300 К, и атомная конфигурация релаксировала к твердому аморфному состоянию (с большим сопротивлением) в процессе молекулярной динамики. Этот основной эффект наблюдался как дальнейший рост числа атомов sp3, и в установившемся стационарном состоянии количество таких атомов составило 10% от общего числа.
5. Электронная структура атомных конфигураций. Проведем анализ электронной структуры эволюции в процессе фазового перехода. Мы будем вычислять ширину запрещенной зоны (gap) Ед(Т) в зависимости от температуры. Для этого мы вычислим электронный спектр системы и определим Ед(Т) = £ьимо ~~ ^номо? гДе HOMO — наивысшая занятая орбиталь, LUMO — низшая незаполненная возбужденная орбиталь, Е — их энергия. Это может быть использовано в рамках многомасштабной модели на основе молекулярной динамики Эренфеста, сформулированной в [13].
Ед(Т), eV
0,6 -
0,4 -
0,2 -
0 1000 2000 3000 4000 Т, К
Рис. 4. Зависимость Ед{Т) для конфигураций, полученных при различных температурах
Зависимость Ед(Т) для конфигураций, полученных при различных температурах, показана на рис. 4. Зависимость ширины энергетической щели от температуры имеет явный немонотонный характер. Первый минимум расположен при комнатной температуре и составляет 0.12-0.16 eV. Второй минимум порядка 0.018 eV в десять раз меньше, чем первый. Эта минимальная ширина щели соответствует температуре Т = 800 К. Поднятие температуры после этой точки приводит к увеличению фракции sp2 атомов. Плато на профиле в диапазоне температур Т = 2000-3000 К соответствует переходу к конфигурации графитовой формы с явно выраженными графитовыми слоями. Около Т = 3700 К графитоподобная конфигурация может быть переведена в графитовую форму с помощью процедуры отжига. Третий минимум порядка 0.08 достигается при температуре Т = 5000 К. Это метастабильная графитоподобная жидкость. После мгновенного охлаждения до температуры Т = 300 К эта конфигурация переходит в алмазоподобную аморфную структуру. Самое большое значение ширины щели 0.68 eV достигается при температуре Т = 600 К. Атомная структура здесь имеет аморфное состояние со связанными графитоподобными и алмазоподобными областями.
6. Заключение. В работе проведено суперкомпьютерное моделирование из первых принципов структурного фазового перехода в аморфном углероде. Моделирование было направлено на объяснение механизма смены сопротивления в экспериментах IBM [3].
Найдены два различных типа атомной структуры. Первый тип структуры, обладающей высоким сопротивлением, появляется при температурах ниже Т = 2000 К. Второй тип конфигурации, обладающей малым сопротивлением, возникает при температурах больше Т = 2000 К. Следует заметить, что изменения в кластерной структуре а-С появляются ниже точки плавления аморфного углерода. Подобные изменения в геометрии системы можно рассматривать как фазовый переход в наноструктуре материала, приводящий к пороговой смене сопротивления.
Авторы благодарят сотрудников Исследовательской лаборатории IBM в Цюрихе докторов Э. Элеф-териу и А. Куриони за инициализацию исследований и постоянное обсуждение результатов вычислений.
СПИСОК ЛИТЕРАТУРЫ
1. Standley В., Bao W., Zhang Н. et al. Graphene-based atomic-scale switches // Nano Lett. 2008. 8. P. 3345-3349.
2. Meijer G. Who wins the nonvolatile memory race? // Science. 2008. 319. P. 1625-1626.
3. Sebastian A., Pauza A., Rossel C., Shelby R.M., Rodriguez A.F., Pozidis H., Elefthe-riou E. Resistance switching at the nanometre scale in amorphous carbon // New J. Phys. 2011. 13. 013020.
4. Иванов-Омский В.И., Лодыгин А.Б., Ястребов С.Г. Сканирующая туннельная спектроскопия аморфного углерода: модель туннелирования // Письма в ЖТФ. 1999. 25. № 24. С. 66-71.
5. Takai К., Oga М., Sato Н. et al. Structure and electronic properties of a nongraphitic disorded carbon system and its heat-treatment effects // Phys. Rev. B. 2003. 67. 214202.
6. Marks N. A., McKenzie D.R., Pailthorpe B.A., Bernasconi M., Parrinello M. Ab initio simulations of tetrahedral amorphous carbon // Phys. Rev. B. 1996. 54. N 14. P. 9703-9714.
7. Yu He, Jinyu Zhang, Ximeng Guan etal. Molecular dynamics study of the switching mechanism of carbon-based resistive memory // IEEE Transactions on Electron Devices. 2010. 57. N 12. P. 3434-3441.
8. Andreoni W., Curioni A. New advances in chemistry and materials science with CPMD and parallel computing // Parallel Computing. 2000. 26. P. 819-842.
9. http://www.cpmd.org.2011
10. Car R., Parrinello M. Unified approach for molecular dynamics and density-functional theory // Phys. Rev. Lett. 1985. 55. P. 2471-2474.
11. Kohn W. Density functional and density matrix method scaling linearly with the number of atoms // Phys. Rev. Lett. 1996. 76. N 17. P. 3168-3171.
12. Bekas C., Curioni A., Andreoni W. Atomic wavefunction initialization in ab initio molecular dynamics using distributed lanczos // Parallel Computing. 2008. 34. N 6. P. 441-450.
13. Shumkin G.N., Popov A.M., Curioni A., Laino T. A multiscale modellling of naphthalocyanine-based molecular switch // Procedia Computer Science. 2010. 1. N 1. P. 185-192.
Поступила в редакцию 07.09.11
52
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2012. № 2
SIMULATION OF RESISTANCE SWITCHING IN AMORPHOUS CARBON
Popov A. M.. Shumkin G. N.
Mathematical simulation of resistance switching mechanism in amorphous carbon (a-C) on the base of quantum molecular dynamics is presented. Electrical conductivity of thin amorphous carbon films is connected to phase transition, induced by electrical field and associated Joule heating and annealing. Simulation shows that transition mechanism is connected to clusterization of existing conductive graphite domains (sp2) inside insulating diamond (sp3) matrix. Conditions of phase transition from one state to another are found from the simulation. Also excited molecular orbitals and energetic gap were calculated at different temperatures.
Keywords: amorphous carbon, phase change material, quantum molecular dynamics, resistance switching, conductivity switching.
