автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Квантовое моделирование переключений в молекулярной структуре наносистем

Московский государственный университет имени М.В. Ломоносова Факультет вычислительной математики и кибернетики

На правах рукописи

Шумкин Георгий Николаевич

Квантовое моделирование переключений в молекулярной структуре наносистем

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

005048104 1 и ЯНВ 2013

Москва - 2012

005048104

Работа выполнена на кафедре автоматизации научных исследований Факультета ВМК МГУ имени М.В. Ломоносова.

Научный руководитель: доктор физико-математических наук,

профессор Попов Александр Михайлович Официальные оппоненты: доктор физико-математических наук,

профессор Мухин Сергей Иванович доктор физико-математических наук, профессор Якобовский Михаил Владимирович

Ведущая организация: Национальный . исследовательский

центр "Курчатовский институт"

Защита состоится IL- 2012 г. вчасов на заседании дис-

сертационного совета Д 501.001.43 при Московском государственном университете имени М.В. Ломоносова, расположенном по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ, 2-й учебный курс, факультет вычислительной математики и кибернетики, ауд. 685

С диссертацией можно ознакомиться в библиотеке факультета ВМК МГУ имени М.В. Ломоносова.

Автореферат разослан

«^L» lÂfatô/lJ 2012 г.

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

диссертационного совета Д 501.001.43, доктор физико-математических наук,

профессор .^—¿bj'L't^Î^ Захаров Е.В.

Общая характеристика работы

Актуальность работы. В последнее время математическое моделирование молекулярных наносиетем приобретает все большее значение, обусловленное необходимостью внедрения новых технологий в микроэлектронике. Современные кремниевые технологии в микроэлектронике имеют целый ряд принципиальных физических ограничений, которые не позволяют более повышать производительность логических схем и увеличивать плотность и быстродействие ячеек памяти. Одним из перспективных подходов в преодолении существующих физических барьеров является развитие молекулярной электроники. Первое направление в молекулярной электронике связано с созданием памяти с сверхвысокой плотностью на основе использования одиночной молекулы в качестве базисной ячейки памяти [1, 2]. Предлагается также использовать одиночную молекулу в качестве диодов, транзисторов и других логических элементов [1,2]. Вторым важным направлением в молекулярной электронике является создание памяти на фазовых переходах (РСМ, Phase Change Memory) [3]. В основе работы такой памяти лежит использование структурного фазового перехода, который происходит под внешним воздействием. Разработка памяти "РСМ" позволит получить значительное повышение плотности хранения информации и быстродействия с сохранением свойства энергонезависимости. Два направления объединяет то, что предлагается создание материалов на основе наносистем, обладающих свойством переключения электропроводности под внешним воздействием (электрического поля, тока и т.д.).

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

описывается квантово-классическими уравнениями движения. Такое моделирование из "первых принципов" стало возможным только в последнее время благодаря использованию современных суперкомпьютеров. Расчеты, выполненные в диссертации, проведены на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ имени М.В. Ломоносова.

Целями диссертационной работы являлись:

1. Разработка многомасштабной квантовомеханической модели для описания переключения электрической проводимости молекулярной структуры.

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

3. Многомасштабное моделирование молекулярного переключателя на основе молекулы нафталоцианина и сравнение с экспериментом.

4. Моделирование структурного фазового перехода в тонкой пленке аморфного углерода на основе термического воздействия на систему. Проведение анализа электронной структуры.

5. Проведение моделирования квантовой молекулярной динамики на суперкомпьютере IBM BlueGene/P. Разработка параллельных алгоритмов и программ для многомасштабной молекулярной динамики.

Научная новизна

1. Сформулирован численный алгоритм вычисления профиля свободной энергии вдоль введенной координационной переменной стохастическим немарковским методом "Метадинамики"

2. Разработана многомасштабная модель переключения электрической проводимости в молекулярной системе. Многомасштабная модель использует вычисленный профиль поверхности свободной энергии, найденный в приближении Борна-Оппенгеймера, и быструю электронную эволюцию системы через возбужденные состояния.

3. Рассчитан профиль поверхности свободной энергии для реакции изомеризации в одиночной молекуле.

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

5. Проведено численное исследование структурного фазового перехода в аморфном углероде и получены энергетические характеристики структур при разных температурах. Проведен анализ температурной эволюции электронной структуры.

мо-

Теоретическая и практическая значимость. Математическое делирование молекулярных систем становится необходимым условием для развития молекулярной электроники, и потому имеет принципиально важное и актуальное значение. Ценность такого квантовомеханического моделирования молекулярной системы из "первых принципов" во многом определяется тем, что эффекты перехода ионов в новую структуру обусловлен взаимодействием с электронной системой, что невозможно моделировать методами классической молекулярной динамики с нейтральными частицами. В диссертации проводится моделирование двух серий экспериментов по созданию прототипов молекулярных ячеек памяти, проведенных в научно-исследовательской лаборатории IBM в Цюрихе [1, 3]. Это определяет практическую

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

Апробация работы. По результатам диссертации были сделаны доклады на следующих конференциях и семинарах:

1. всероссийской суперкомпыотерной конференции "Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность" (г. Новороссийск, 21-26 сентября 2009 г);

2. XVII-ой международной научной конференции студентов, аспирантов и молодых ученых "Ломоносов-2010", секция "Вычислительная математика и кибернетика" (г. Москва, 12-15 апреля 2010);

3. международной конференции по вычислительным наукам ICCS 2010 (International Conference on Computational Science 2010, Нидерланды, г. Амстердам, 31 мая - 2 июня 2010);

4. международной конференции по вычислительным наукам ICCS 2012 (International Conference on Computational Science 2012, США, шт. Небраска, г. Омаха, 4-6 июня 2012);

5. семинаре "Advancing the Frontiers of Modeling and Simulation in Chemistry and Materials Science, IBM Zurich Laboratory" (Швейцария, г. Цюрих, январь 2011).

Публикации. Основные результаты диссертации опубликованы в одиннадцати работах: из них шесть статей в журналах, рекомендованных ВАК, [Al, А2, A3, А4, А5, А6], три статьи в тематических сборниках [А7, А8, А9], два тезиса докладов [А10, All]

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

Содержание работы

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

Первая глава посвящена математической постановке квантовой молекулярной динамики ("основная задача") в приближении Борна-Оппенгейме-ра. Приводится разработанная в диссертации многомасштабная модель. Постановке "основной задачи" и постановке многомасштабной модели посвящены параграфы § 1.1 и § 1.2, соответственно.

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

= (1)

где — координаты ионов, М/ — массы составляющих их ядер, N1 — число ионов 1-го типа, Гц, — силы, действующие на ионы со стороны электронов и других ионов.

Расчет сил на ионы для основного квантового состояния осуществляется на основе вычисления электронной плотности р(г) как функции трех пере-

менных и координат ионов И/:

(¿Г,

(2)

II/ - г

где Zl — заряды ионов.

Для нахождения электронной плотности р(т) необходимо найти орби-тали ^п(г) Кон-Шэма [4], являющиеся решением следующей стационарной задачи для уравнения Шредингера на собственные значения при заданном положении ионов

где ие//(р, г) — эффективный потенциал, зависящий от электронной плотности р и положения ионов И/, еп — энергия п-ой орбитали.

Эффективный потенциал [/е//(р, г) представляет собой сумму

где [1ХС — обменно-корреляционный потенциал, представляющий собой аналитическую зависимость от электронной плотности р и градиента плотности Vр в приближении РВЕ [5].

Электростатический потенциал Хартри [6] £/# определяется из уравнения Пуассона

Потенциал 1]ц представляет собой коллективный потенциал электронов, тогда как обменный-корреляционный член 11хс содержит все многочастичные электронные эффекты.

Электронная плотность представляется в виде разложения по одноэлек-тронным орбиталям Кон-Шэма фп:

(3)

Дин = -47гр(г).

где Ne — количество электронов.

Нелинейность задачи на собственные значения связана с тем, что эффективный потенциал Ucff(p,r) зависит от электронной плотности р(г), которая в свою очередь зависит от решения задачи — электронных орбиталей фп. Поэтому необходим итерационный процесс. Допустим, что на начальной итерации известно некоторое приближение для плотности. По плотности определяется пространственная зависимость эффективного потенциала Ueff(p, г). Тогда используя решение задачи на собственные значения еп и собственные функции -фп, вычисляется новая пространственная зависимость электронной плотности системы по формуле (6). Сходимость процесса определяется нормой разницы между новой и старой плотностями.

В диссертации для решения нелинейной задачи на собственные значения на каждом шаге квантовой молекулярной динамики используются метод Кар-Парринелло [7] и созданный на его основе код "CPMD" (Car-Parrinello Molecular Dynam¡cs)[8], Уравнения квантовой молекулярной динамики Кар-Парринелло выглядят следующим образом:

Qlф /j, 2 '

М-= ^ ~ Ueff(P> г' + Anmí'm( Г, t)

т—1

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

■фп(г, t)ipm(r, t)dr = Snm, (8)

где 5nm — символ Кронекера. В уравнении для электронов в (7) вводится фиктивная масса электрона fi, которая определяет скорость релаксации электронной системы к основному состоянию при каждом новом положении ионов. Множители кпт представляют собой химический потенциал, возникающий

9

из условия сохранения количества электронов в системе, что приводит к требованию ортонормированности (8) орбиталей Кон-Шэма {Ф„}. Преимущество метода Кар-Парринелло заключается в том, что можно один раз решить задачу на собственные значения, а далее основное состояние находится с помощью искусственного итерационного процесса. Итерационный процесс получается из использования явной центрированной схемы при дискретизациии временного уравнения для электронных орбиталей фп в (7).

Задача квантовой молекулярной динамики решается для многих ионов, и цель определения динамики частиц состоит в том, чтобы после многих шагов по времени усреднением по времени и по числу частиц определить термодинамические параметры и их соотношения между собой - уравнения состояния системы.

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

Множество А состоит из группы соседних ионов, чье поведение определяет

(9)

где относительное расстояние гу определяется согласно

(10)

процесс переключения атомной структуры.

Для моделирования свойства переключения в системе в диссертации разработана многомасштабная модель, которая позволяет описывать процесс переключения с использованием переходов системы через возбужденные состояния. Для этого в рамках временной теории функционала плотности одноэлек-тронная волновая функция раскладывается по полному набору собственных функций — орбиталей Кон-Шэма г!>п (3):

ф (Г, {К/},0 = (г, {К/}). (11)

п

Построение мпогомасштабной модели для моделирования свойства переключения включает следующие шаги и приближения:

1. Вводится координационная переменная в реакции переключения. Координационная переменная характеризует переход из одной структуры в другую.

2. Нахождение профиля поверхности свободной энергии Гельмгольца вдоль выбранной координационной переменной в стохастическим немарковским методом "Метадинамики" в рамках решения "основной задачи".

3. Вычисление набора необходимых возбужденных орбиталей Кон-Шэма Фп-

4. Решение временной электронной задачи на основе вычисленного профиля свободной энергии.

Постановка задачи для многомасштабной модели с использованием ко-

ординационной переменной принимает следующий вид:

^=- £ - Е £ сыь* - «¿«и«)

п>т0 т>та п>тц

1Ьп = е^К - г ^ Ст«и.(з)а + ^ УтшгСт' ^ т>то т>то

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

Уравнения в (12) связаны через коэффициенты (1пт{з):

= (13)

Коэффициенты К™ определяющие действие внешнего потенциала, вводятся следующим образом:

\ф*п(т)ЩтШг, (14)

где V — оператор внешнего потенциала.

Задача при интегрировании уравнений (12) - определить временную динамику электронных коэффициентов с„ и координационной переменной а.

В параграфе § 1.3 представлены численные методы для решения "основной" задачи и численная схема для интегрирования уравнений многомасштабной модели. Для численного решения нелинейной задачи на собственные функции используется базисное разложение функций фп в ряд Фурье:

где И - объем элементарной ячейки, к - вектор обратной решетки и с71,к -Фурье-коэффициенты орбитали п. При представлении решения в виде ряда Фурье оператор кинетической энергии в (3), представляющий собой сумму вторых производных, приводит просто к умножению на квадрат волнового

12

числа. Умножение орбитали фп на эффективный потенциал проводится в реальном пространстве. Для преобразования величин между двумя пространствами используется быстрое преобразование Фурье.

Для интегрирования уравнения Кар-Парринелло [7] используется численная схема Верлета, для ортогонализации орбиталей в методе Кар-Парринелло на каждом шаге интегрирования используется алгоритм SHAKE/RATTLE.

Для решения уравнения Пуассона (5) при периодических граничных условий используется быстрое преобразование Фурье. Для граничных условий, соответствующих моделированию изолированной молекулы или изолированного кластера, используется сочетание методов функции Грина и быстрого преобразования Фурье.

Для численного интегрирования уравнений (12) многомасштабной модели предлагаются неявная схема для электронного уравнения и явная центрированная схема для уравнения координационной переменной.

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

Параграф § 1.5 посвящен тестированию параллельного кода CPMD. В параграфе представлено исследование параллельных алгоритмов кода CPMD на примере расчетов молекулы нафталоцианина, и найдены оптимальные параметры распараллеливания кода CPMD для использования на суперкомпь-тере IBM BlueGene/P. Показано, что параллельная эффективность кода может достигать 70%.

В параграфе § 1.6 представлены исследования численных параметров молекулярной динамики Кар-Парринелло. Проанализировано, сколько членов в ряде Фурье требуется для проведения расчетов электронной структуры с достаточной точностью, для чего были проведены вычисления потенциальной энергии системы и ширины запрещенной зоны. Определена граница устойчивости численной схемы Кар-Парринелло.

Во второй главе представлены результаты численного моделирования молекулярного переключателя на основе одиночной молекулы нафталоциа-нина. Геометрия молекулы изображена на рис. 1(а). Механизм переключения в молекуле нафталоцианина основывается на реакции водородной изомеризации, в результате которой два атома водорода во внутренней полости молекулы (рис. 1(а)) меняют свое положение, что приводит к образованию второго изомера (рис. 1(6)). Реакция изомеризации в эксперименте запускается с помощью тока сканирующего туннельного микроскопа. Такая реакция равносильна повороту всей молекулярной структуры на 90°, что ведет в свою очередь к смене электрической проводимости вдоль двух перпендикулярных осей симметрии молекулы нафталоцианина. Особенность представленного молекулярного переключателя состоит в том, что оба изомера имеют одну геометрию, в отличие от ранее исследованных.

В параграфе § 2.1 представлена математическая постановка моделирования химической реакции изомеризации в молекулярном переключателе методом "Метадинамики". Метод "Метадинамики" представляет собой мощный стохастический немарковский метод восстановления поверхности свободной энергии Гельмгольца вдоль выбранной координационной переменной. Метод "Метадинамики" основан на введении в систему искусственного зависящего от времени потенциала ^(в,^, который добавляется в систему в определен-

ные моменты времени tf

u<t

где s1 = s(ti). Этот потенциал V есть сумма гауссианов, каждый из которых имеет высоту Я¿ и полуширину W{. Потенциал V выталкивает частицу из ямы функции F (s). Профиль поверхности свободной энергии F(s) в в конце расчета "Метадинамики" должен быть равен взятому с обратным знаком потенциалу V(s,t):

F(s) =-V(s,t). (17)

Трудность расчетов с применением стохастического метода "Метадинамики" связана с тем, что шаг "Метадинамики" tmtd = U+ï — U не может быть произвольно малым. Фактически, он на порядок и более превосходит временной шаг квантовой молекулярной динамики. Для расчетов необходимо выбрать оптимальные параметры tmtd, щ, щ с точки зрения точности вычислений и времени вычислений. Сложность использования метода "Метадинамики" состоит в нахождении последовательности включения потенциалов, приводящей к сходимости метода.

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

Параграф § 2.3 посвящен вычислениям профиля поверхности свободной энергии F(s) методом "Метадинамики" и определению величины энергетического барьера реакции изомеризации. В параграфе определены численные

ж}

(16)

параметры Я,-, тмто метода "Метадинамики" для достижения сходимости. Проведен анализ временной эволюции координационной переменной я, и найдены три потенциальные ямы, в которых может находиться переменная й. Выбор шага "Метадинамики" тмто связан с необходимостью релаксации системы после добавления каждого Гауссиана. Последовательные расчеты с уменьшающейся высотой Гауссианов Я, при сохранении' остальных параметров "Метадинамики" постоянными позволили установить сходимость при вычислении профиля поверхности свободной энергии и величины энергетического барьера для химической реакции водородной изомеризации в молекуле нафталоцианина. Вычисленные профили поверхности свободной энергии вдоль выбранной координационной переменной в последовательных расчетах "Метадинамики" с уменьшающейся высотой Гауссианов представлены на рис. 2. Найдено, что энергетический барьер реакции водородной изомеризации, полученный методом "Метадинамики" равен 52 кДж/моль.

Параграф § 2.4 посвящен исследованию влияния геометрии молекулярного переключателя на высоту энергетического барьера реакции изомеризации. Найдено, что добавление дополнительных ароматических колец к рукавам молекулы нафталоцианина (рис. 1(а)) позволяет снизить барьер химической реакции изомеризации с 52 кДж/моль до 35 кДж/моль. Это связывается с тем, что увеличение длины молекулы уменьшает связи атомов водорода внутри полости молекулы с атомами азота.

В параграфе § 2.5 проведен анализ эволюции электронной плотности в течение процесса "Метадинамики". Выяснено, что эволюция электронной плотности молекулы нафталоцианина на различных радиусах от центра молекулы существенно определяет характер протекания реакции изомеризации. Этот вывод важен для построении многомасштабной модели, в которой главной (координационной) переменной должна стать угловая переменная при определенном радиусе от центра молекулы.

а

б

в г

Рис. 1. Оптимизированная геометрия молекулы нафталоцианина для двух изомеров и вычисленные возбужденные орбитали для первого изомера: (а) начальная конфигурация молекулы нафталоцианина, (б) конечная конфигурация, обладающая геометрией второго изомера, (в) рассчитанная структура молекулярных орбиталей LUMO и (г) LUMO t1. Белые шары соответствуют ядрам водорода, черные - ядрам азота, серые ядрам углерода. LUMO (Lowest Unoccupied Molecular Orbital) — низшая незанятая молекулярная орбиталь.

о

-60

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 S

Рис. 2. Вычисленный профиль F(s) поверхность свободной энергии вдоль координационной переменной s. Сходимость достигнута в результате четырех последовательных расчетов с параметрами : 1 )Н = ЗкТ; W = 0.03 2) Я = 1.5fcT;3) Я = 0.75fcT;4) Я = 0.375/сГ, где к — константа Вольцмана, Т — температура (ЗОК)

Параграф § 2.6 посвящен расчетам квантовой многомасштабной модели для молекулярного переключателя и сравнению результатов расчетов с экспериментом. В эксперименте наблюдалось сильное различие скоростей переключения при разных положениях иглы сканирующего туннельного микроскопа (СТМ) над молекулой нафталоцианина. Введем следующую координатную систему: молекула нафталоцианина располагается в плоскости XY, координата плоскости соответствует координате zq вдоль оси г. Центр молекулы имеет координаты (xq, уо, го). Потенциал, создаваемый иглой СТМ, моделировался с помощью следующей формы:

где го соответствует плоскости молекулы, Г] и сто — положение и эффективный радиус иглы СТМ, соответственно; Уо определяет величину напряжения. Для точек, находящихся на разном удалении от центра молекулы, найдены

V = У0 • (г - z0)e(r-T>)2/'

(18)

пороговые напряжения, при которых происходит переключение. Для наиболее удаленной точки пороговое напряжение является минимальным среди всех остальных точек. В центральной точке для максимально рассмотренного порогового напряжения переключение не наблюдается. Полученные результаты расчетов многомасштабной модели находятся в хорошем согласии с экспериментальными наблюдениями. Пример модельного расчета для наиболее удаленной от центра молекулы точки приведен на рис. 3, 4. Для расчета используется набор рассчитанных возбужденных орбиталей LUMO и LUMO+l((LUMO (Lowest Unoccupied Molecular Orbital) — низшая незанятая молекулярная орбиталь)), геометрия которых дана на рис. 1(в) и рис. 1(г), соответственно, а также невозбужденная орбиталь HOMO ((Highest Occupied Molecular Orbital) — высшая занятая молекулярная орбиталь). Орбитали LUMO и LUMO-f 1 совпадают с точностью до поворота на 90°. Таким образом, электронное переключение возбужденных орбиталей LUMO и LUMO+1 должно приводить к переключению системы. Перед началом расчета коэффициент для HOMO равен 1, для других орбиталей 0. В начале расчета происходит быстрое переключение между невозбужденной орбиталью HOMO и возбужденными орбиталями LUMO и LUMO+1 (рис. 3). В конечном счете переключение электронных орбиталей HOMO, LUMO и LUMO+1 на рис. 3 приводит к переключению переменной s(t) на рис. 4. Переключение s(t) соответствует перемещению атомов водорода в молекуле нафталоциаиина иод действием внешнего потенциала, создаваемого иглой СТМ, что и есть основа процесса переключения.

В третьей главе представлены результаты моделирования структурного фазового перехода в тонкой пленке аморфного углерода.

Параграф § 3.1 посвящен методике проведения расчетов для моделирования структурного фазового перехода в аморфном углероде. Начальный состав элементарной ячейки строился по следующему принципу: в алмазной

40000

80000

Time, at. un.

Рис. 3. Электронное переключение. HOMO (Highest Occupied Molecular Orbital) — высшая занятая молекулярная орбиталь. LUMO (Lowest Unoccupied Molecular Orbital) — низшая незанятая молекулярная орбиталь.

2.5

1.5

0.5

Ш

— sW

„iijl

««see

1 II

ЩЦ

40000

80000

Time, at. un.

Рис. 4. Соответствующее переключение переменной s. 20

Рис. 5. Геометрия молекулярной системы после процесса молекулярной динамики, проведенной при (а) Г = 2000Л: и (б) Г = 3000К.

матрице располагался графитовый кластер в качестве начальной затравки. Далее проводились расчеты квантовой молекулярной динамики при различных температурах до установления термодинамического равновесия. Затем к полученным атомным конфигурациям применялся отжиг. Расчеты показали, что молекулярную структуру существенно определяют сценарии поднятия температуры и отжига. Показано, что при повышении температуры и последующем отжиге получаются проводящие графитовые структуры. При температурах порядка Т = 300(Ж в окончательной структуре не остается алмазоподобньтх атомов, что видно на рис. 5, где сравниваются геометрии, полученные при 200(Ж и 300(Ж. Кристаллическая структура при Т = ЗОООЛГ получается независимо от начальных условий. При дальнейшем повышении температуры до 5000К и последующем быстром отжиге получается аморфное состояние из кристаллического. Это соответствует тому, что в эксперименте наблюдается как переход из непроводящего аморфного состояния в проводящее кристаллическое, так и обратный переход.

В параграфе § 3.2 для оценки смены проводимости приведен анализ электронной структуры атомных конфигураций аморфного углерода, полу-

Рис. 6. Температурная зависимость ширины запрещенной зоны Ед{Т) для аморфного углерода.

ченных при разных температурах. Для оценки электропроводности найдена температурная зависимость ширины запрещенной зоны для аморфного углерода, представленная на рис. 6. Т.к. атомные конфигурации в расчетах имеют аморфную структуру, то используется статистический подход для подсчета ширины запрещенной зоны. Дисперсия ширины запрещенной зоны Ед(Т) показана на рис. 6 вертикальными отрезками "Error". Также для оценки электропроводности проведены расчеты при различных температурах уровня Ферми и плотности энергетических состояний на уровне Ферми.

Параграф § 3.3 посвящен результатам по построению перколяционных структур. Разработана программа на основе алгоритмов перколяции для анализа электронной плотности как функции трех переменных. Найдена геометрия проводящих графитовых кластеров на основе анализа электронной плотности атомных конфигураций, полученных в молекулярной динамике при различных температурах.

Заключение. В диссертации предложены модели квантовой молекуляр-

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

1. Разработана математическая модель молекулярного переключателя в одиночной молекуле нафталоцианина на основе молекулярной динамики Борна-Оппенгеймера и Эренфеста. Модель позволяет рассчитывать профиль свободной энергии вдоль выбранной координационной переменной и переход системы в новую проводящую фазу через возбужденные энергетические состояния.

2. Разработаны параллельные алгоритмы и программы, реализующие расчеты многомасштабной молекулярной динамики на суперкомпьютере IBM BlueGene/P.

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

4. Проведено исследование процесса структурного фазового перехода в аморфном углероде. Найдено, что перестройка молекулярной струкгу-

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

Благодарность Автор выражает благодарность своему научному руководителю, доктору физико-математических наук, профессору Александру Михайловичу Попову за постоянное внимание к работе.

Список публикаций

[А1] Г.Н Шумкин, A.M. Попов, А. Куриони, Т. Лайно. Моделирование из первых принципов молекулярного переключателя на основе реакции изомеризации // Мателшпическое моделирование, — 2010.— Vol. 22, по. 11.-Pp. 18-28.

[А2] G.N. Shumkin, A.M. Popov, A. Curioni, T. Laino. A multiscale modelling of naphthalocyanine-based molecular switch // Procedia Computer Science. - 2010. - Vol. 1, no. 1. - Pp. 185-192.

[A3] A.M. Попов, Г.Н. Шумкин. Многомасштабная модель молекулярного переключателя // Вестн. Моск. Ун-та, сер. 15, Вычисл. матем. и киберн. - 2010. - по. 3. - Pp. 11-18.

[А4] Г.Н. Шумкин, A.M. Попов. Моделирование из первых принципов фазового перехода в аморфном углероде // Математическое моделирование. — 2012. — Vol. 24, по. 10. — Pp. 65-79.

[А5] G.N. Shumkin, F. Zipoli, A.M. Popov, A. Curioni. Multiscale Quantum

Simulation of Resistance Switching in Amorphous Carbon // Procedía Computer Science. — 2012. — Vol. 9. — Pp. 641-650.

[A6] A.M. Попов, Г.Н. Шумкин. Моделирование переключения сопротивления в аморфном углероде // Вестпн. Моск. Ун-та, сер. 15, Вычисл. матем. и киберп. — 2012. — по. 2. — Pp. 3-9.

[А7] Г.Н. Шумкин. Численное ( моделирование молекулярных нанопереключателей // Тезисы лучших дипломных работ факультета ВМК МГУ им. Ломоносова. - 2009. — Р. 39.

[А8] Г.Н. Шумкин. Математическое моделирование молекулярных переключателей // "Программные системы и инструменты": Тематический сборник факультета ВМК МГУ им. Ломоносова. -2009,- Vol. 10.- Pp. 140-149.

[А9] Г.Н. Шумкин. Суперкомпьютерное моделирование памяти на основе фазового перехода // "Программные системы и инструменты": Тематический сборник факультет,а ВМК МГУ им. Ломоносова. -2011,-Vol. 12.-Pp. 99-110.

[AlO] H.H. Попова, O.B. Докосан, Г.Н. Шумкин. Методы визуальной поддержки для задач молекулярного моделирования на суперкомпьютере BlueGene/P // Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность: Труды Всероссийской суперкомпьютерной конференции (21-26 сентября 2009 г., г. Новороссийск). — 2009. — Pp. 426-428.

[All] Г.Н. Шумкин. Вычисления из первых принципов пути реакции изомеризации в молекулярной системе методом метадинамики // Сборник тезисов XVII Международной научной конференции

студентов, аспирантов и молодых ученых Ломоносов-2010, секция Вычислительная математика и кибернетика. — 2010. — Р. 152.

Цитированная литература

[1] P. Liljeroth, J. Repp, G. Meyer. Current-Induced Hydrogen Tautomeriza-tion and Conductance Switching of Naphthalocyanine Molecules // Science. — 2007. - Vol. 317. - Pp. 1203-1206.

[2] N.J. Tao. Electron transport in molecular junctions /,/ Nanotech.— 2006. — Vol. l.-Pp. 173-181.

[3] A. Sebastian, A. Pauza, C. Rossel et al. Resistance switching at the nanometre scale in amorphous carbon /,/ New J. Phys. - 2011. - Vol. 13. - P. 013020.

[4] W. Kohn, L. J. Sham. Self-Consistent Equations Including Exchange and Correlation Effects // Phys. Rev. - 1965. - Vol. 140, no. 4A. - P. A1133-A1138.

[5] J. P. Perdew, K. Burke, M. Ernzerhof. Generalized Gradient Approximation Made Simple // Phys. Rev. Lett. - 1996. - Vol. 77, no. 18. - Pp. 3865-3868.

[6] D. Marx, J. Hutter. Ab initio molecular dynamics: Theory and implementation 11 Modern Methods and Algorithms of Quantum Chemistry. — 2000.— Vol. l.-Pp. 329-477.

[7] R. Car, M. Parrinello. Unified approach for molecular dynamics and density-functional theory ,// Phys. Rev. Lett. - 1985. - Vol. 55,- Pp. 2471-2474.

[8] The CPMD Consortium,http://www.cpmd.org, 2012.

Напечатано с готового оригинал-макета

Подписано в печать 22.11.2012 г. Формат 60x90 1/16. Усл.печ.л. 1,0. Тираж 100 экз. Заказ 458.

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к. Тел. 8-495-939-3890. Тел./факс 8-495-939-3891.

Введение

Глава 1. Постановка задачи.

1.1. Математические модели "из первых принципов".

1.2. Учет быстрой релаксации электронов в переключении структуры атомов. Многомасштабная модель переключения

1.3. Численные методы решения задачи

1.4. Параллельные алгоритмы решения основной задачи

1.5. Тестирование параллельного кода СРМБ.

1.6. Исследование свойств численного метода.

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

2.1. Математическая постановка моделирования реакции изомеризации

2.2. Определение координационной переменной.

2.3. Определение профиля поверхности свободной энергии

2.4. Способы понижения барьера.

2.5. Анализ электронной плотности

2.6. Многомасштабная модель переключателя

Глава 3. Моделирование переключения в нанопленке аморфного углерода

3.1. Вычисление энергетической кривой для фазового перехода

3.2. Анализ электронной структуры

3.3. Построение перколяционных структур

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

В последнее время математическое моделирование молекулярных наносистем приобретает все большее значение, обусловленное необходимостью внедрения новых технологий в микроэлектронике. Современные кремниевые технологии в микроэлектронике имеют целый ряд принципиальных физических ограничений, которые не позволяют более повышать производительность логических схем и увеличивать плотность и быстродействие ячеек памяти. Одним из перспективных направлений в преодолении существующих физических барьеров является развитие молекулярной электроники. Для создания энергонезависимой памяти с чрезвычайно высокой плотностью предлагается использовать одиночную молекулу в качестве базисной ячейк памяти [1-5]. Кроме того, предлагается использовать одиночную молекулу в качестве диодов, транзисторов и других логических элементов [1, 6-9]. Также важным направлением в молекулярной электронике является также создание памяти на фазовых переходах(РСМ, Phase Change Memory), в основе которой лежат тонкие пленки халькогенидных стекол. Толщина пленок составляет несколько десятков нанометров. В основе работы такой памяти лежит использование структурного фазового перехода [10, 11], который происходит под внешним воздействием. Также исследовались различные углеродные материалы (углеродные нанотрубки, графен) в качестве материала для памяти на фазовых переходах [12]. Развитие этого направления также обещает получить значительное повышение плотности хранения информации и быстродействия с сохранением свойства энергонезависимости. Эти направления объединяет то, что в обоих случаях предлагается создание материалов на основе наносистем, обладающих свойством переключения электропроводности под внешним воздействием (электрического поля, тока и т.д.). Переключение электропроводности в сотни и тысячи раз происходит в областях, имеющих размер порядка 1-100 нм.

Диссертация посвящена численному квантовомеханическому моделированию молекулярных систем с переключательным поведением. Молекулярная система представляется набором ионов и электронов, совместная эволюция которых описывается квантово-классическими уравнениями движения. Такое моделирование из "первых принципов" стало возможным только в последнее время благодаря использованию современных суперкомпьютеров. Расчеты, выполненные в диссертации проведены на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ имени М.В. Ломоносова. Ценность такого "прямого" моделирования молекулярной системы во многом определяется сравнением результатов моделирования с экспериментом.

В диссертации рассматриваются две физические проблемы, связанные с исследованием подобных наносистем.

Первая проблема связана с исследованием молекулярного переключателя на основе одиночной молекулы нафталоцианина методами математического моделирования. Молекулярным переключателем называют молекулярное соединение, которое обладает двумя или более устойчивыми состояниями, отличающимися электропроводностью на несколько порядков. Переключение между состояниями представляет собой прохождение молекулярной системы через энергетический барьер поверхности свободной энергии. Устойчивые состояния соответствуют минимумам поверхности свободной энергии. В общем случае переключение, т.е. прохождение через барьер, может быть вызвано действием тока с использованием установки сканирующего туннельного микроскопа (СТМ), действием электрического поля, действием света, а также возможно термическое возбуждение системы, вызывающее переход в новое состояние. Возможно также использование окислительно-восстановительных химических реакций для изменения проводимости молекулярной системы, для чего используется установка электрохимического СТМ. Исчерпывающее описание механизмов работы различных типов молекулярных переключателей дано в обзоре [13]. В этой же работе описаны основные проблемы использования молекулярных переключателей в индустрии.

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

Исследование молекулярного переключателя на основе реакции изомеризации в молекулярной системе мотивировано последними успешными экспериментами, выполненными в цюрихской исследовательской лаборатории 1ВМ[14]. Переключение проводимости в молекуле нафталоцианина осуществляется на основе реакции водородной изомеризации (таутомеризации). Во внутренней полости молекулы два атома водорода могут менять свое положение в молекуле, что приводит к изменению проводимости через молекулу под действием электрического потенциала сканирующего туннельного микроскопа. При этом перемещение двух атомов водорода не меняет геометрию молекулы, а приводит к переходу от одного изомера к другому. Такая реакция равносильна повороту всей молекулярной структуры на 90°, что ведет в свою очередь к смене молекулярной проводимости вдоль двух перпендикулярных осей симметрии молекулы нафталоцианина. Эксперименты показали, что процесс изомеризации проходит без разрушения геометрии молекул в отличие от большинства ранее предложенных молекулярных переключателей, основанных на разрушении геометрии молекулярных систем при переключении [2, 15-17].

Химическая реакция, приводящая к смене проводимости, запускается в результате действия тока Сканирующего Туннельного Микроскопа "СТМ". Над мономолекулярной тонкой пленкой нафталоцианина, напыленной на соляную подложку, под которой находится медная подложка, располагается игла СТМ. При подаче достаточно большого напряжения (1.7В) в установке СТМ сила тока начинает переключаться между двумя хорошо определёнными уровнями. Подробное исследование, представленное в работе [14], установило однозначную связь между переключениями тока, т.е. переключениями молекулярной проводимости, с реакцией изомеризации. Более того, было установлено, что наибольшая скорость переключения достигается, когда игла располагается на периферии молекулы, более чем в 10 А от центра молекулы. Также было найдено, что вероятность переключения заметно различалась на двух различных рукавах молекулы (подробное описание строения молекулы будет дано в Главе 2). Представленный в экспериментах [14] молекулярный переключатель хорошо подходит для использования технологии самособирающихся монослоев, что позволяет говорить о масштабируемости подобной молекулярной технологии. Однако прямому выходу экспериментов в технологическую стадию мешает ряд проблем, связанных с наблюдаемыми эффектами. Так, при подаче определенного потенциала над соответствующей точкой молекулы молекулярная система начинает переключаться и в соседних полостях. Таким образом, наблюдается нелокальный эффект переключения. Также для молекулярного переключателя необходимо определить, при какой энергии возникает химическая реакция таутомеризации. Кроме того, требуется объяснить экспериментальную зависимость скорости переключения от положения иглы над молекулой нафталоциапииа и ответить на вопрос, почему скорость переключения на периферии молекулы наивысшая.

На выяснение причин такого поведения и направлены исследования, проводимые в диссертации.

Вторая рассматриваемая в диссертации физическая проблема связана с численным моделированием структурного фазового перехода в тонких пленках аморфного углерода, приводящего с существенному переключению проводимости системы. Эксперименты направлены на создание памяти на фазовых переходах на основе использования аморфного углерода. Интерес к этой задаче возник в связи с недавними успешными экспериментами, проведенными в цюрихской исследовательской лаборатории IBM [18]. Исследователями предложен аморфный углерод в качестве нового материала для памяти на фазовых переходах. Ранее в качестве подобного рассматривались в основном халь-когеиидные стекла, которые исследовались на протяжении длительного времени. В экспериментах [18] обнаружено и исследовано переключение проводимости в тонких пленках аморфного углерода. С помощью атомносилового микроскопа (AFM, Atomic Force Microscope) исследователям удалось получить переключение электрической проводимости аморфного углерода в нано-точке, диаметром 20 нм, в тонких пленках с толщиной несколько десятков нанометров. Переключение проводимости происходило при определенном пороговом напряжении. Величина порогового напряжения зависела линейно от толщины пленки и составляла 2В для толщины пленки 20 нм. При регистрации тока во время этого процесса переключения получалась ¿"-образная вольт-амперная характеристика. Подобные 5-образные вольт-амперные характеристики наблюдаются при переключении проводимости в тонких пленках на основе халькогенидных стекол [10, 11]. Это экспериментальное наблюдение позволяет рассматривать аморфный углерод в качестве перспективного материала для РСМ. Во-первых, это связано с тем, что технология приготовления аморфного углерода хорошо известна и сравнительно дешева. Во-вторых, открытие нового материала для памяти на фазовых переходах дает возможность значительно расширить спектр исследований в этой области.

В работе [18] экспериментаторами были сформулированы теоретические проблемы, которые необходимо разрешить до внедрения экспериментов в индустрию. Требуется ответить на вопрос, каков механизм, ответственный за процесс переключения. Экспериментаторы обнаружили, что 5-образные вольт-амперные характеристики хорошо приближаются кривыми, получаемые с помощью модели Пула-Френкеля[19], однако используемое значение диэлектрической константы отличается в сто раз от физических значений для аморфного углерода. Выдвигаемые исследователями гипотезы механизма переключения основаны на температурных, полевых (электрических) эффектах, а также кластеризации графитоподобных атомов углерода. Также экспериментаторами поставлен вопрос, может ли регибридизация атомов углерода быть ответственной за механизм переключения электропроводности в тонкой пленке аморфного углерода.

Методы исследования

Выбор методов математического моделирования связан напрямую с вопросами, которые поставлены в экспериментах [14, 18]. Из-за того, что в основе вышеизложенных экспериментальных наблюдений лежат квантовые эффекты, требуется использовать квантовомеханические методы, которые в соединении с молекулярной динамикой ионов позволяют вычислять поверхность свободной энергии Гельмгольца динамической системы, а значит, вычислять барьер химической реакции изомеризации в молекулярном переключателе наиболее точным способом, позволяют учитывать рассчитывать впрямую электронную гибридизацию атомов углерода, исследовать влияние внешнего электрического поля и температуры. Используемые в диссертации модели называются моделями из "первых принципов". В отличие от полуэмпирических методов, где используются эмпирические (подгоночные) параметры в потенциале взаимодействия между частицами, в вычислениях из "первых принципов" включено только кулоновское взаимодействие.

Использование квантовомеханических методов прямого расчета молекулярной системы фиксированной конфигурации ядер основывалось в начале на методе Хартри-Фока [20], сформулированном уже в начале ЗОх годов 20-го века. Применение этих методов было сильно ограничено по различным причинам. Во-первых, для достижения т.н. "химической точности "в случае системы без симметрии необходимо минимально рассматривать порядка ТУо ^ Ю(с точностью до порядка множителя порядка 2) атомов [21]. Для методов, основанных на многочастичных волновых функциях, таких, как метод Хартри-Фока, даже на современных вычислительных системах эта величина является пределом. Сложность расчета взаимодействия частиц в методе Хартри-Фока масштабируется как ./V3 — И4. При моделировании молекулярных систем со свойством переключения электропроводности потребовалось бы моделировать порядка нескольких сотен атомов для получения правильных статистик в квантовой молекулярной динамике для различных конфигураций ядер. Квантовое моделирование систем, состоящих из сотен атомов стало возможно, только с момента формулировки теории функционала плотности [22, 23] на суперкомпьютерах. За разработку теории функционала плотности и за появление первых кодов квантовой химии присуждена Нобелевская премия по химии В. Кону и Дж. Поплу в 1998 г. Теория функционала плотности позволяет свести многоэлектронную квантовую (Ы-мерную) задачу для основного состояния к поиску трехмерной электронной плотности. Во-вторых, начавшийся бурный рост суперкомпыотерной отрасли в середине 80х годов 20-го века позволил начать проводить квантовое моделирование на основе теории функционала плотности. Постоянное увеличение производительности производимых суперкомпьютеров, идущее и по сегодняшний день, позволило на текущий момент моделировать квантовые наносистемы, состоящие из сотен и даже тысяч атомов [21, 24-27]. Все расчеты, результаты которых приведены в диссертации, проведены на суперкомпьютере IBM BlueGene/P. Один вариант расчета на этом суперкомпьютере в диссертационных исследованиях занимал порядка 20 часов, что эквивалентно 3 годам расчета на самом быстром однопроцессорном компьютере. Соединение теории функционала плотности и методов молекулярной динамики позволило существенно продвинуться в расчетах молекулярных систем. Последующее продвижение также связано с появлением метода Кара-Парринелло [28]. В методе введена искусственная динамика электронов, позволяющая не проводить расчет стационарной задачи на собственные значения на каждом шаге перемещения ионов. Электронная система релаксирует к новому положению ионов за счет искусственной электронной динамики, введенной в систему. Молекулярная динамика Кара-Парринелло оказалась с вычислительной точки зрения очень успешной и потому нашла применение в различных физических, химических и биологических задачах [29]. На сегодняшний день созданы различные стандартные программы квантовой молекулярной динамики[24, 25, 30] на основе теории функционала плотности. Во всех этих программах также используется теория псевдопотенциалов [29, 31, 32], в которой поле остовных электронов и ядра заменяется псевдопотенциалом. Это приближение связано с тем, что в химических реакциях участвуют только валентные электроны. Для расчетов использовался квантовомеханический код молекулярной динамики CPMD (Car-Parrinello Molecular Dynamics, Молекулярная динамика Кара-Парринелло) [24], разработанный в Цюрихской исследовательской лаборатории. Создание такого кода заняло по крайней мере два десятка лет, код был протестирован на множестве молекулярных систем. Многолетнее тестирование кода позволило ведущим журналам мира принимать статьи, содержащие результаты расчетов с использованием кода СРМБ [24]. Однако несмотря на то, что сегодня присутствуют подобные параллельные стандартные программы, позволяющие проводить на суперкомпьютерах квантовоме-ханические расчеты для сотен и тысяч атомов, при их создании использовались, конечно, ряд сильных приближений. Прежде всего, это то, что системы могут моделироваться только для основного квантового состояния. В диссертации рассматривались наносистемы, обладающие свойством переключения электропроводности. Для исследования таких наносистем в диссертации разработана многомасштабная модель, которая позволяет учитывать переход системы из одной транспортной фазы в другую через возбужденные состояния. Из расчетов квантовой молекулярной динамики с помощью кода СРМБ вычисляется профиль свободной энергии вдоль координационной переменной. Координационная переменная описывает главное направление процесса переключения атомной структуры, приводящего к смене электропроводности, и ее введение позволяет свести вычисление многомерной поверхности свободной энергии в общем случае к вычислению одномерного профиля поверхности свободной энергии. Более строгое определение координационной переменной представлено в главе 1. С помощью кода СРМБ также вычисляются возбужденные орбитали, перекрытие которых определяет процесс переключения атомной структуры через возбужденные состояния. Для построения многомасштабной модели используется вычисленный профиль поверхности свободной энергии и вычисленные возбужденные орбитали. Описание уравнений многомасштабной модели также представлено в главе 1.

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

В коде СРМБ для интегрирования уравнений Кар-Парринелло используется метод Фурье: электронные орбитали раскладываются по обратным векторам элементарной ячейки (плоским волнам). При представлении решения в виде ряда Фурье оператор кинетической энергии в одноэлектронном уравнении Шредингера, представляющий собой сумму вторых производных, приводит просто к умножению на квадрат волнового числа. Умножение электронной орбитали на эффективный потенциал проводится в реальном пространстве. Для преобразования величин между двумя пространствами используется быстрое преобразование Фурье, распараллеливание которого определяет эффективность кода СРМЭ. В основе распараллеливания лежит равномерное распределение плоских волн Фурье-пространства по процессорам и двухмерных плоскостей в реальном пространстве.

Цели исследования

Цели исследования состояли в следующем:.

1. Разработка многомасштабной квантовомеханической модели для описания переключения электрической проводимости молекулярной структуры.

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

3. Многомасштабное моделирование молекулярного переключателя на основе молекулы нафталоцианина и сравнение с экспериментом.

4. Моделирование структурного фазового перехода в тонкой пленке аморфного углерода на основе термического воздействия на систему. Проведение анализа электронной структуры.

5. Проведение моделирования квантовой молекулярной динамики на суперкомпьютере IBM BlueGene/P. Разработка параллельных алгоритмов и программ для многомасштабной молекулярной динамики.

Научная новизна

1. Сформулирован численный алгоритм вычисления профиля свободной энергии вдоль введенной координационной переменной стохастическим немарковским методом "Метадинамики"

2. Разработана многомасштабная модель переключения электрической проводимости в молекулярной системе. Многомасштабная модель использует вычисленный профиль поверхности свободной энергии, найденный в приближении Борна-Оппенгеймера, и быструю электронную эволюцию системы через возбужденные состояния.

3. Рассчитан профиль поверхности свободной энергии для реакции изомеризации в одиночной молекуле.

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

5. Проведено численное исследование структурного фазового перехода в аморфном углероде и получены энергетические характеристики структур при разных температурах. Проведен анализ температурной эволюции электронной структуры.

Практическая и теоретическая значимость результатов.

Математическое моделирование молекулярных систем становится необходимым условием для развития молекулярной электроники, и потому имеет принципиально важное и актуальное значение. Расчеты, выполненные в диссертации, проведены на суперкомпьютере IBM Blue Gene/P, установленном на факультете ВМК МГУ имени М.В. Ломоносова. Ценность такого "прямого" квантовомеханического моделирования молекулярной системы во многом определяется тем, что эффекты перехода ионов в новую структуру обусловлен взаимодействием с электронной системой, что невозможно моделировать методами классической молекулярной динамики с нейтральными частицами.

Тем не менее полезность разрабатываемых моделей определяется прежде всего сравнением результатов моделирования с экспериментом. В диссертации проводится моделирование двух серий экспериментов по созданию молекулярных переключателей, проведенных в научно-исследовательской лаборатории IBM в Цюрихе [14, 18]. Это определяет практическую ценность моделирования, т.к. указанные эксперименты связаны с выходом на новые технологии в области создания сверхплотной памяти с высоким быстродействием. Практическая важность связана также с тем, что моделирование переключения в молекулярных системах носит принципиально квантовомеханическое. Использование эмпирических потенциалов взаимодействия в молекулярной динамике не позволит моделировать процесс переключения между двумя состояниями, имеющих различную атомную структуру. Эмпирический потенциал пригоден только для описания одной равновесной структуры. В кван-товомеханической модели не используются подгоночные параметры, поэтому она способна описывать любую атомную структуру из "первых принципов а также переходы между ними, т.е. процесс переключения. Квантовомеханиче-ская модель учитывает также не только надбарьерные переходы системы, но и эффекты туннелирования сквозь потенциальный барьер.

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

Апробация работы

По результатам диссертации были сделаны доклады на следующих конференциях:

1. всероссийской суперкомпьютерной конференции "Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность" (г. Новороссийск, 21-26 сентября 2009 г);

2. XVII-ой международной научной конференции студентов, аспирантов и молодых ученых "Ломоносов-2010", секция "Вычислительная математика и кибернетика" (г. Москва, 12-15 апреля 2010);

3. международной конференции по вычислительным наукам ICCS 2010(1п-ternational Conference on Computational Science 2010, Нидерланды, г. Амстердам, 31 мая - 2 июня 2010);

4. международной конференции по вычислительным наукам ICCS 2012(1п-ternational Conference on Computational Science 2012, США, шт. Небраска, г. Омаха, 4-6 июня 2012);

5. семинаре "Advancing the Frontiers of Modeling and Simulation in Chemistry and Materials Science, IBM Zurich Laboratory" (Швейцария, г. Цюрих, январь 2011).

Публикации

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

1. Г.Н. Шумкин, A.M. Попов, А. Куриони, Т. Лайно. Моделирование из первых принципов молекулярного переключателя на основе реакции изомеризации // Математическое моделирование. // 2010. Т. 22, № 11. С. 18-28.

2. G.N. Shumkin, A.M. Popov, A. Curioni, Т. Laino. A multiscale modelling of naphthalocyanine-based molecular switch // Procedia Computer Science. //2010. Vol. 1, no. 1. Pp. 185-192.

3. A.M. Попов, Г.Н. Шумкин. Многомасштабная модель молекулярного переключателя // Вестн. Моск. Ун-та, сер. 15, Вычисл. матем. и киберн. 2010. № 3. С. 11-18.

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

5. G.N. Shumkin, F. Zipoli, A.M. Popov, A. Curioni. Multiscale Quantum Simulation of Resistance Switching in Amorphous Carbon // Procedia Computer Science.// 2012. Vol. 9. Pp. 641-650.

6. A.M. Попов, Г.Н. Шумкин. Моделирование переключения сопротивления в аморфном углероде // Вестн. Моск. Ун-та, сер. 15, Вычисл. матем. и киберн. 2012. № 2. С. 3-9.

7. Г.Н. Шумкин. Численное моделирование молекулярных нанопереклю-чателей // Тезисы лучших дипломных работ факультета ВМК МГУ им. Ломоносова. 2009. С. 39.

8. Г.Н. Шумкин. Математическое моделирование молекулярных переключателей // "Программные системы и инструменты": Тематический сборник факультета ВМК МГУ им. Ломоносова. 2009. Т. 10. С. 140-149.

9. Г.Н. Шумкин. Суперкомпьютерное моделирование памяти на основе фазового перехода // "Программные системы и инструменты": Тематический сборник факультета ВМК МГУ им. Ломоносова. // 2011. Т. 12. С. 7-14.

10. H.H. Попова, О.В. Джосан, Г.Н. Шумкин. Методы визуальной поддержки для задач молекулярного моделирования на суперкомпьютере BlueGene/I // Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность: Труды Всероссийской суперкомпьютерной конференции (21-26 сентября 2009 г., г. Новороссийск). 2009. С. 426-428.

11. Г.Н. Шумкин. Вычисления из первых принципов пути реакции изомеризации в молекулярной системе методом метадинамики // Сборник тезисов XVII Международной научной конференции студентов, аспирантов и молодых ученых Ломоносов-2010, секция Вычислительная математика и кибернетика. 2010. С. 152.

Структура диссертации

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

Заключение

Заключение. В диссертации предложены модели квантовой молекулярной динамики для молекулярных систем, обладающих переключательным свойством проводимости, основанным на изменении молекулярной структуры. Разработана программа для многомасштабного моделирования переключательных свойств. Проведено моделирование двух процессов переключения в структуре молекулярной системы: процесса водородной изомеризации в одиночной молекуле пафталоцианина и структурного фазового перехода в аморфном углероде. Расчеты проведены на суперкомпьютере IBM BlueGene/P. Получены следующие основные результаты:

1. Разработана математическая модель молекулярного переключателя в одиночной молекуле нафталоцианина на основе молекулярной динамики Борна-Оппенгеймера и Эренфеста. Модель позволяет рассчитывать профиль свободной энергии вдоль выбранной координационной переменной и переход системы в новую проводящую фазу через возбужденные энергетические состояния.

2. Разработаны параллельные алгоритмы и программы, реализующие расчеты многомасштабной молекулярной динамики на суперкомпьютере IBM BlueGene/P.

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

4. Проведено исследование процесса структурного фазового перехода в аморфном углероде. Найдено, что перестройка молекулярной структуры из аморфной в кристаллическую связана с термическими эффектами. Получены энергетические характеристики системы в различных структурных состояниях. Показано изменение электропроводности при фазовом переходе на основе эволюции электронной структуры и плотности квантовых состояний системы на уровне Ферми.

1. Tao, N. Electron transport in molecular junctions / N. Tao // Nanotech. — 2006.-Vol. 1.- Pp. 173-181.

2. A 2]catenane-based solid state electronically reconfigurable switch / C. Collier, J. Stoddart, J. Heath, et. al. // Science2000.- Vol. 289.-Pp. 1172-1175.

3. Controlling the charge state of individual gold adatoms / J. Repp, G. Meyer, F. Olsson, M. Persson // Science. 2004. - Vol. 305,- Pp. 493-495.

4. Szuchmacher Blum, A. Molecularly inherent voltage-controlled conductance switching / A. Szuchmacher Blum, et. al. // Nature Materials. — 2005. — Vol. 4. Pp. 167-172.

5. Rieth, M. Handbook of Theoretical and Computational Nanotechnology / M. Rieth, W. Schommers. — Karlsruhe, Germany: Forschungszentrum Karlsruhe, 2006.

6. Eigler, D. M. An atomic switch realized with the scanning tunnelling microscope / D. M. Eigler, C. P. Lutz, W. E. Rudge // Nature1991.- Vol. 352. Pp. 600-603.

7. Joachim, C. Electronics using hybrid-molecular and mono-molecular devices / C. Joachim, J. K. Gimzewski, A. Aviram // Nature. — 2000.— Vol. 408. Pp. 541-548.

8. Nitzan, A. Electron transport in molecular wire junctions / A. Nitzan, M. A. Ratner // Science. 2003. - Vol. 300. - Pp. 1384-1389.

9. Selzer, Y. Single-molecule electrical junctions / Y. Selzer, D. Allara // Annu. Rev. Phys. Chem. 2006. - Vol. 57. - Pp. 593-623.

10. Meijer, G. Who wins the nonvolatile memory race? / G. Meijer // Science. — 2008. Vol. 319. - Pp. 1625-1626.

11. Wuttig, M. Phase-change materials for rewritable data storage / M. Wuttig, N. Yamada // Nat. Mater. 2007. - Vol. 6. - Pp. 824-832.

12. Kreupl, F. Carbon-based resistive memory / F. Kreupl, R. Bruchhaus, et. al. // IEEE Int. Electron Devices Meeting Technical Digest. — 2008.— Pp. 521-524.

13. Liljeroth, P. Charge transport through molecular switches / P. Liljeroth, S. Molen // Journal of Physics: Condensed Matter. — 2010.— Vol. 22, no. 13,- P. 133001.

14. Liljeroth, P. Current-Induced Hydrogen Tautomerization and Conductance Switching of Naphthalocyanine Molecules / P. Liljeroth, J. Repp, G. Meyer // Science. 2007. - Vol. 317. - Pp. 1203-1206.

15. Conformational changes of single molecules induced by scanning tunneling microscopy manipulation: A route to molecular switching / F. Moresco, G. Meyer, K. Rieder et al. // Phys. Rev. Lett. 2001. - Vol. 86. - Pp. 672-675.

16. Qiu, X. H. Mechanisms of reversible conformational transitions in a single molecule / X. H. Qiu, G. V. Nazin, W. Ho // Phys. Rev. Lett. 2004.-Vol. 93. - Pp. 196806-196809.

17. Iancu, V. Realization of a four-step molecular switch in scanning tunneling microscope manipulation of single chlorophyll-a molecules / V. Iancu, H. S.W. // Proc. Natl. Acad. Sci. U.S.A.- 2006,- Vol. 103.-Pp. 13718-13721.

18. Resistance switching at the nanometre scale in amorphous carbon / A. Sebastian, A. Pauza, C. Rossel et al. // New J. Phys. — 2011.— Vol. 13.— P. 013020.

19. Frenkel, J. On Pre-Breakdown Phenomena in Insulators and Electronic Semi-Conductors / J. Frenkel // Phys. Rev.- 1938,- Vol. 54, no. 8.— Pp. 647-648.

20. A.M. Попов. Вычислительные нанотехнологии / A.M. Попов. — Москва: Макс Пресс, 2009.

21. Кон, В. Электронная структура вещества — волновые функции и функционалы плотности. (Нобелевская лекция. Стокгольм. 28 января 1999 года.) / В. Кон // Успехи Физических Наук, — 2002.— Vol. 172, по. 3. Pp. 336-348.

22. Kohn, W. Density Functional and Density Matrix Method Scaling Linearly with the Number of Atoms / W. Kohn // Phys. Rev. Lett. 1996. - Vol. 76, no. 17.- Pp. 3168-3171.

23. Kohn, W. Self-Consistent Equations Including Exchange and Correlation Effects / W. Kohn, L. J. Sham // Phys. Rev. 1965.- Vol. 140, no. 4A.-P. A1133-A1138.

24. The CPMD Consortium,http://www.cpmd.org, 2012.

25. The Regents of the University of California,http://eslab.ucdavis.edu/, 2012.

26. Andreoni, W. New Advances in Chemistry and Materials Science with CPMD and Parallel Computing / W. Andreoni, A. Curioni // Parallel Computing. — 2000. Vol. 26. - Pp. 819-842.

27. Gygi, F. Large-scale first-principles molecular dynamics: moving from teras-cale to petascale computing / F. Gygi // Journal of Physics: Conference Series. 2006. - Vol. 46. - Pp. 268-277.

28. Car, R. Unified approach for molecular dynamics and density-functional theory / R. Car, M. Parrinello 11 Phys. Rev. Lett.- 1985.- Vol. 55.-Pp. 2471-2474. ^i

29. Marx, D. Ab initio molecular dynamics: Theory and implementation / D. Marx, J. Hutter // Modern Methods and Algorithms of Quantum Chemistry. 2000. - Vol. 1. - Pp. 329-477.

30. The Vienna Ab initio Simulation Package, http://www.vasp.at/, 2012.

31. Topp, W. Chemically Motivated Pseudopotential for Sodium / W. Topp, J. John J. Hopfield // Phys. Rev. B. 1973. - Vol. 7, no. 4. - Pp. 1295-1303.

32. Th. Starkloff, T. Local pseudopotential theory for transition metals / T. Th. Starkloff, J. Joannopoulos // Phys. Rev. B.- 1977,- Vol. 16, no. 12,- Pp. 5212-5215.

33. Stm-induced switching of the hydrogen molecule in naphthalocyanine / S. A., N. Arboleda, M. David, et. al. // J. Phys.: Condens. Matter. — 2009.— Vol. 21.- P. 064201.

34. Fu, Q. Mechanism for tautomerization induced conductance switching of naphthalocyanine molecule / Q. Fu, J. Yang, Y. Luo // Appl. Phys. Lett. — 2009,-Vol. 95,- P. 182103.

35. Henkelman, G. A climbing image nudged elastic band method for finding saddle points and minimum energy paths / G. Henkelman, B. Bias P. Uberuaga, H. Jonsson // J. Chem. Phys. 2000. - Vol. 113, no. 22,- Pp. 9901-9904.

36. Kumar, S. The weighted histogram analysis method for free-energy calculations on biomolecules. i: the method / S. Kumar, D. Bouzida, R. Swendsen // J. Comput. Chem. 1992. - Vol. 13, no. 8.- Pp. 1011-1021.

37. Kastner, J. Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: Umbrella integration / J. Kastner, W. Thiel // J.Chem. Phys.- 2005,- Vol. 123, no. 14.-P. 144104.

38. A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactions / B. Ensing, A. Laio, M. Parrinello, et. al. // J. Phys. Chem. B. 2005. - Vol. 109, no. 14. - Pp. 6676-6687.

39. Robertson, J. Diamond-like amorphous carbon / J. Robertson // Materials Science and Engineering: R: Reports. — 2002. — Vol. 37. — Pp. 129-281.

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

41. Thermal stability and rehybridization of carbon bonding in tetrahedral amorphous carbon / D. Grierson, A. Sumant, A. Konicek, et. al. // J. Appl. Phys. 2010. - Vol. 107. - P. 033523.

42. Stress reduction and bond stability during thermal annealing of tetrahedral amorphous carbon / A. Ferrari, B. Kleinsorge, M. N.A., et. al. // J. Appl. Phys. 1999. - Vol. 85. - Pp. 7191-7197.

43. Иванов-Омский, В. Сканирующая туннельная спектроскопия аморфного углерода: модель туннелирования / В. Иванов-Омский, А. Лодыгин, С. Ястребов // Письма в ЖТФ. 1999. - Vol. 25, по. 24,- Pp. 66-71.

44. Ab initio simulations of tetrahedral amorphous carbon / N. A. Marks, D. R. McKenzie, B. A. Pailthorpe et al. // Phys. Rev. B. 1996. - Vol. 54, no. 14. - Pp. 9703-9714.

45. Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon / J. Tersoff // Phys. Rev. Lett. 1988.- Vol. 61, no. 25,- Pp. 2879-2882.

46. TersoffJ. Modeling solid-state chemistry: Interatomic potentials for multi-component systems / J. Tersoff // Phys. Rev. B. — 1989. — Vol. 39, no. 8. — Pp. 5566-5568.

47. Molecular dynamics study of the switching mechanism of carbon-based resistive memory / Y. He, J. Zhang, X. Guan, et. al // IEEE Transactions on Electron Devices. 2010. - Vol. 57, no. 12. - Pp. 3434-3441.

48. Binder, K. Monte Carlo Simulation in Statistical Physics / K. Binder,

49. D. W. Heermann. — Springer, 2010.

50. Perdew, J. P. Generalized Gradient Approximation Made Simple / J. P. Perdew, K. Burke, M. Ernzerhof // Phys. Rev. Lett. 1996. - Vol. 77, no. 18. - Pp. 3865-3868.

51. Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions / W. G. Hoover // Phys. Rev. A. 1985. - Vol. 31, no. 3.- Pp. 1695-1697.

52. Martyna, G. J. Nose-Hoover chains: The canonical ensemble via continuous dynamics / G. J. Martyna, K. M. L., M. Tuckerman // J. Chem. Phys.— 1992. Vol. 97, no. 4. - Pp. 2635-2643.

53. Runge, E. Density-functional theory for time-dependent systems / E. Runge,

54. E. Gross // Phys. Rev. Lett. — 1984. — Vol. 52, no. 12.- Pp. 997-1000.

55. Verlet, L. Computer "Experiments"on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules / L. Verlet // Phys. Rev. — 1967. — Vol. 159, no. 1. Pp. 98-103.

56. Three real-space discretization techniques in electronic structure calculations. / T. Torsti, T. Eirola, J. Enkovaara, et.al // physica status solidi (b). 2006. - Vol. 243, no. 5. - Pp. 1016-1053.

57. Finite-element implementation for electron transport in nanostructures / P. Havu, V. Havu, M. Puska, et. al. // J. Chem. Phys.- 2006.- Vol. 124. P. 054707.

58. Bekas, C. Atomic wavefunction initialization in ab initio molecular dynamics using distributed lanczos / C. Bekas, A. Curioni, W. Andreoni // Parallel Computing. 2008. - Vol. 34, no. 6. - Pp. 441-450.

59. Fukui, K. The path of chemical reactions, the IRC approach / K. Fukui // Acc. Chem. Res. 1981. - Vol. 14. - Pp. 363-368.

60. Lonfat, M. The energy gap of carbon clusters studied by scanning tunneling spectroscopy / M. Lonfat, B. Marsen, К. Sattler // Chem. Phys. Lett.— 1999. Vol. 313. - Pp. 539-543.

61. Kubo, R. Statistical mechanics: an advanced course with problems and solutions / R. Kubo. North Holland, 1999.

62. Kohn, W. Inhomogeneous Electron Gas / W. Kohn, P. Hohenberg // Phys. Rev. 1966. - Vol. 136. - Pp. B864-B871.

63. Хокни, P. Численное моделирование методом частиц / Хокни, Р., Иствуд, Дж. — Москва: Мир, 1987. — Р. 640.