Математическое моделирование характеристик автоэлектронной эмиссии наноструктур

Сорокина Вероника Андреевна

Математическое моделирование характеристик автоэлектронной эмиссии наноструктур

Изучение эмиссионных характеристик наноструктур чрезвычайно важно в перспективе создания новых типов электронных приборов. Предполагается, что использование полевых эмиттеров на основе нанотрубок позволит получать значительные в сравнении с другими источниками эмиссионные токи при относительно невысоких значениях приложенного напряжения. Такая возможность существует благодаря высокому соотношению их геометрических размеров. Работа посвящена теоретическому исследованию влияния внешнего электрического поля на различные идеальные структуры открытых одностенных нанотрубок углерода и карбида кремния. Методом решения является теория функционала электронной плотности с градиентным приближением для обменно-корреляционного функционала (B3LYP). Начальное приближение сформировано в базисе 6-31g. Влияние внешнего электрического поля на объекты изучения оценивалось по изменениям распределения электронной плотности в структурах, плотности тока эмиссии, полных электронных энергий, дипольных моментов и ширины запрещённых зон. Из теоретических расчётов были найдены значения тока эмиссии, соответствующие каждой модели. Качественный характер изменения плотности эмиссионного тока отвечает ожидаемому по экспериментальным данным результату. С точки зрения эмиссионных характеристик наиболее выгодной оказалась нанотрубка карбида кремния. Также на основе анализа зонной структуры объектов сделан вывод, что вне зависимости от конфигурации углеродные нанотрубки должны обладать полупроводниковым типом проводимости.

Общественные науки в целом

Диссертации

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36825f1be77c40d5910d

UUID: 248c5a5f-03ce-44a1-b7d6-697f901fd639

Язык: Русский

Опубликовано: больше 7 лет назад

Просмотры: 30

Сорокина Вероника Андреевна

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 1417050 bytes

Поделиться работой

Санкт-Петербургский государственный университет Факультет прикладной математики — процессов управления Кафедра моделирования электромеханических и компьютерных систем Сорокина Вероника Андреевна Магистерская диссертация Математическое моделирование характеристик автоэлектронной эмиссии наноструктур Направление 03.04.01 Прикладные математика и физика Магистерская программа «Прикладная информатика» Научный руководитель, кандидат физ.-мат. наук, доцент Никифоров К. А. Санкт-Петербург 2016

Содержание 1. Введение 4 2. Постановка задачи 6 3. Предмет и объекты изучения 9 3.1. Полевая электронная эмиссия . . . . . . . . . . . . . . . . . . 9 3.2. Нанотрубки . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2.1. Углеродные нанотрубки . . . . . . . . . . . . . . . . . 12 3.2.2. Карбидокремниевые нанотрубки . . . . . . . . . . . . 16 4. Методы исследования 19 4.1. Математическая модель задачи . . . . . . . . . . . . . . . . . 19 4.2. Обзор методов решения уравнения Шрёдингера . . . . . . . 22 4.3. Математический аппарат . . . . . . . . . . . . . . . . . . . . 24 4.3.1. Вариационный принцип квантовой механики . . . . . 24 4.3.2. Вариационный метод Ритца . . . . . . . . . . . . . . . 25 4.4. Особенности теоретических методов исследования . . . . . . 27 4.5. Вычислительные методы квантовой химии . . . . . . . . . . 28 4.5.1. Метод Хартри . . . . . . . . . . . . . . . . . . . . . . . 28 4.5.2. Метод Хартри—Фока . . . . . . . . . . . . . . . . . . . 29 4.5.3. Метод функционала электронной плотности . . . . . 32 5. Практическая реализация 43 5.1. Формирование волновой функции . . . . . . . . . . . . . . . 43 5.1.1. Плоские волны . . . . . . . . . . . . . . . . . . . . . . 44 5.1.2. Орбитали Слэтера . . . . . . . . . . . . . . . . . . . . 45 5.1.3. Орбитали Гаусса . . . . . . . . . . . . . . . . . . . . . 47 5.2. Численное интегрирование выражений для вычисления обменно-корреляционной энергии . . . . . . 50 5.3. Запуск задач в Gaussian и распределённые вычисления . . . 52 2

6. Результаты численного эксперимента 56 6.1. Результаты моделирования структур . . . . . . . . . . . . . . 56 6.2. Распределение электронной плотности . . . . . . . . . . . . . 58 6.3. Плотность эмиссионного тока . . . . . . . . . . . . . . . . . . 61 6.4. Полные энергии и дипольные моменты систем . . . . . . . . 63 6.5. Анализ зонной структуры нанотрубок . . . . . . . . . . . . . 64 7. Заключение 67 3

1. Введение Нанотрубки могут быть рассмотрены не только как источники зондирующих электронных пучков, когерентных в случае полевой электронной эмиссии, но и как объекты зондирования этими пучками [1]. В обоих случаях наличие внешнего электрического поля обуславливает необходимость моделирования влияния этого поля на характеристики нанотрубок. Анализ существующих исследований в области полевой электронной эмиссии показывает, что материал автоэмиссионного катода, работающего в высоком техническом вакууме, должен обладать специфической комбинацией свойств, а именно, низким и постоянным значением работы выхода электронов в сочетании с высокой механической прочностью и долговечностью, а также высокой тепло- и электропроводностью. Кроме того, также требуется, чтобы эти материалы были технологически доступны. В связи с этим, наиболее потенциально перспективными семействами материалов являются материалы на основе углерода и карбида кремния, так как именно они показывают преимущества во всех перечисленных выше свойствах, помимо требований доступности. Это подтверждают текущие данные по полевой эмиссии нано- [2] и мезомасштабных [3] структур карбида кремния и углерода [4]. Синтез углеродных и карбидокремниевых нанотрубок с заданными свойствами в количестве, достаточном для массового производства, в настоящее время не налажен, что обостряет актуальность теоретического моделирования поведения данных структур в сильном электрическом поле. Апробация работы По материалам работы представлены доклады в рамках XLV, XLVI, XLVII международных научных конференций аспирантов и студентов «Процессы управления и устойчивость» (Control Processes and Stability — CPS’14, CPS’15, CPS’16), III международной конференции «Устойчивость и процессы управления» (Stability and Control Processes in Memory of V. I. Zubov — SCP’15), V международной конференции «Cовременные тенденции научных исследований нанообъектов искусственного и природного 4

происхождения» (State-of-the-art Trends of Scientific Research of Artificial and Natural Nanoobjects — STRANN’16). Структура и объём работы Работа состоит из введения, пяти основных разделов и заключения. Во втором разделе сформулирована цель работы, поставлена задача исследования и оговорены применяемые методы решения. В третьей части проведен обзор литературы, позволяющий установить актуальность темы исследования. В четвертой главе излагаются теоретические основы, на которых базируется работа квантово-химических программ. Пятый раздел посвящен вопросам практической реализации задач, в шестом представлены оригинальные результаты исследования. В заключении сформулированы достижения и выводы. Работа изложена на 71 листах печатного текста, содержит 2 блоксхемы, 25 рисунков, 2 графика и 6 таблиц. Библиографический список состоит из 36 ссылок. 5

2. Постановка задачи Целью работы является теоретическое исследование влияния внешнего электрического поля на различные структуры открытых одностенных нанотрубок углерода и карбида кремния. Величина напряжённости поля выбрана исходя из характерного для явления полевой электронной эмиссии значения. В качестве моделей выбраны две углеродные нанотрубки с одинаковым количеством атомов (обе прямые), но различных диаметров — C72 (3,0) (рис. 2.1) и C72 (3,3) (рис. 2.2); карбидокремниевая нанотрубка Si36 C36 (3,3) (рис. 2.3) той же длины, и две углеродные нанотрубки одной конфигурации, но разных длин: C120 (5,5) (рис. 2.4) и C240 (5,5) (рис. 2.5). Выбор объектов изучения обусловлен стремлением определить наиболее подходящую конфигурацию и длину наноструктур с точки зрения эмиссионных характеристик. Рис. 2.1: Углеродная нанотрубка C72 (3, 0) Рис. 2.2: Углеродная нанотрубка C72 (3, 3) Эти материалы имеют чётко организованную структуру. В качестве метода моделирования выбрана теория функционала электронной плотности. С её применением связаны вопросы способа задания начального приближения и обменно-корреляционного функционала. Средой моделирования структур является программный пакет для квантово-химических 6

Рис. 2.3: Карбидокремниевая нанотрубка Si36 C36 (3, 3) Рис. 2.4: Углеродная нанотрубка C120 (5, 5) Рис. 2.5: Углеродная нанотрубка C240 (5, 5) вычислений Gaussian версии 09. Вычисления были проведены на кластере высокопроизводительных вычислений факультета прикладной математики — процессов управления и на виртуальной машине ресурсного центра Научного парка СПбГУ «Вычислительный центр». Для достижения указанной цели необходимо: • провести обзор литературы по углеродным и карбидокремниевым нанотрубкам, определить их место в развитии новых типов электронных приборов, отметить преимущества перед другими источниками и сравнить между собой; • обосновать применение метода функционала электронной плотности, указать альтернативные методы исследования и отметить значение 7

выбранного метода среди них; • привести возможные способы задания пробной волновой функции и обменно-корреляционного функционала, выбрать подходящее приближение; • изучить схему взаимодействия с программой Gaussian 09 и кластером высокопроизводительных вычислений; • получить координатное представление необходимых структур и «поместить» их в сильное внешнее электрическое поле; • сравнить результаты моделирования структур вне и под воздействием поля, оценить эмиссионные характеристики нанотрубок различных конфигураций. 8

3. Предмет и объекты изучения Широта возможных применений материала, обещающего стать революционным, породила огромное количество исследований. Нанотрубкам прочат большое будущее и в микроэлектронике (диоды, полевые транзисторы), и в оптике (дисплеи, светодиоды), и в механике (сверхпрочные нити, композитные материалы), и в микроскопии (зонды сканирующих туннельных и атомно-силовых микроскопов), и в медицине (безопасная транспортировка медикаментов и осуществление новой технологии лечения заболеваний) [5]. Уникальные электрические и механические свойства нанотрубок делают их привлекательными и для использования в качестве полевых эмиттеров. Полевые эмиттеры на основе нанотрубок обладают узким энергетическим спектром, что особенно важно для создания дисплеев с высоким разрешением и микроскопии, а также помогают избежать больших затрат энергии (например, на нагрев) и потенциального ущерба материалу, быстро реагируя на изменение внешнего электрического поля [6]. По сравнению с другими источниками эмиттеры из нанотрубок позволяют получать более высокие эмиссионные токи при относительно невысоких значениях приложенного напряжения. Причиной тому служит их необычная атомная структура. 3.1. Полевая электронная эмиссия Электронной эмиссией называется испускание электронов некоторой поверхностью. Электроны в материале ограничены потенциалом, и для того, чтобы покинуть металл, им необходимо преодолеть потенциальный барьер. Это может произойти в результате нагрева (термоэлектронная эмиссия), облучения (фотоэлектронная эмиссия), бомбардировки электронами, ионами (вторичная электронная и ионно-электронная эмиссии соответственно), агрессивного воздействия (экзоэлектронная эмиссия) или в случае нали9

чия у поверхности тела сильного электрического поля. Последнее явление именуется автоэлектронной или полевой электронной эмиссией. Её главное преимущество перед другими видами состоит в отсутствии предварительного возбуждения электронов: для реализации не требуются дополнительные затраты энергии [7]. Электроны приобретают способность выходить в вакуум только за счёт искривления формы потенциального барьера. Это явление впервые было описано в 1928 г. физиками Р. Фаулером и Л. Нордгеймом, хотя наблюдалось ещё в конце XIX столетия. В то время оно не могло быть объяснено в силу отсутствия аппарата квантовой механики. В исследовании плотности тока эмиссии особенно ценна традиционная модель полевой эмиссии — теория Фаулера—Нордгейма, позволяющая получить значение эмиссионного тока на основании напряжённости приложенного поля. Результат подтверждается многими экспериментами, однако некоторые наноэмиттеры составляют исключение [8]. Это связано с особенностью нанотехнологий. Их целью является создание такой структуры материала на атомном уровне, которая существенно изменила бы характерные физические процессы на пространственных масштабах. Таким образом становится возможным получение улучшенных характеристик системы, которое, в то же время, и усложняет описание процесса [9]. В теории Фаулера—Нордгейма предполагается, что электроны не взаимодействуют друг с другом, но для наноэмиттеров межэлектронным взаимодействием пренебречь нельзя. Эта квантово-механическая черта полевой эмиссии выходит далеко за рамки теории Фаулера—Нордгейма и не вполне изучена экспериментально или теоретически. Теория функционала электронной плотности в применении к наноэмиттерам позволяет решить эту проблему. Следует учитывать, что во время автоэлектронной эмиссии катод разогревается из-за разницы между средней энергией электронов, подходящих к поверхности катода, и средней энергией электронов, уходящих сквозь потенциальный барьер (т. н. эффект Ноттингема). Кроме того, в процессе эмиссии существует вероятность повреждения эмиттера вслед10

ствие действия того же сильного поля. Это большая проблема, которая, тем не менее, может быть решена выбором специального материала. Также полевую эмиссию можно наблюдать при меньших напряжениях, если придать катоду форму тонкого острия. На сегодняшний день наиболее устойчивыми к таким повреждениям считаются углеродные нанотрубки [6]. Изучение аллотропных модификаций углерода вообще представляет большой интерес, так как позволяет перейти на преимущественно новый уровень использования материалов, прочных на растяжение и изгиб. 3.2. Нанотрубки Исторической предпосылкой начала изучения нанотрубок явилось открытие четвертой аллотропной формы углерода (помимо графита, алмаза и карбина; см. рис. 3.1) — фуллерена — нобелевскими лауреатами Р. Кёрлом, Х. Крото и Р. Смолли в 1985 г. Фуллерен представляет собой шарообразное молекулярное соединение, состоящее из чётного числа атомов углерода. Рис. 3.1: Аллотропные модификации углерода 11

Хотя упоминания и предсказания о существовании образований, схожих с нанотрубками, появлялись ранее, первые были описаны лишь в 1991 г. одновременно двумя группами учёных: японскими под руководством С. Ииджимы и русскими во главе с Л. А. Чернозатонским. Материал был открыт при исследовании продуктов, образующихся при разряде вольтовой дуги в атмосфере гелия. Речь идёт о кристаллах, атомы которых формируют полый стержень, диаметр варьируется от одного до нескольких десятков нанометров, а длина достигает нескольких микрон. Первоначально возможность образования частиц в виде трубок была обнаружена для углерода, затем для нитрида бора, карбида кремния и других соединений [6]. На сегодняшний день известно об устойчивости нанотрубок к критическому механическому воздействию и деформациям, вызванным изменениями температуры или радиоактивным излучением [10]. Они способны к самовосстановлению: при возникновении повреждения ячейки циклически передвигаются вдоль поверхности трубки, перераспределяя энергию. Кроме того, они демонстрируют множество неожиданных электрических, магнитных, оптических свойств, которые уже стали объектами ряда исследований [6], [5], [11]. Особенностью углеродных нанотрубок является высокая по сравнению с другими известными проводниками электропроводность, теплопроводность и химическая стабильность. 3.2.1. Углеродные нанотрубки Углеродная нанотрубка представляет собой свёртку графенового листа (рис. 3.2). Графен — двухмерная плоская молекула углерода толщиной в один атом (рис. 3.1) — был получен позже углеродных нанотрубок, вследствие чего изучен слабее. Повышенная мобильность электронов переводит его в разряд наиболее перспективных материалов для наноэлектроники [12]. Свёртка графенового листа в однородный по длине цилиндр может происходить вдоль различных направлений. Поэтому однослойные трубки характеризуются хиральным вектором (n, m), где числа n и m (индексы хиральности) являются координатами вектора в кристаллографиче12

Рис. 3.2: Условная схема свёртки графенового листа ском базисе графенового листа и однозначно определяют диаметр нанотрубок (рис. 3.3). Через индексы хиральности диаметр нанотрубки можно найти по формуле √ d= 3d0 p 2 m + n2 + mn, π где d0 = 1,42 Å — расстояние между соседними атомами углерода в графеновой плоскости. На практике измеряется именно эта величина, так как разрешающая способность современных электронных микроскопов не позволяет различать трубки разных типов симметрий. Долгое время считалось, что индексы хиральности определяют также и зонную структуру: если разность n и m кратна трём или n = m, то материал ведёт себя как проводник, иначе — как полупроводник. Такое мнение, вероятно, устоялось после выхода в 1992 г. статей группы американских учёных, получивших одними из первых энергетический спектр углеродных нанотрубок [13]. Однако в последнее время стали появляться работы [5], ставящие этот факт под сомнение. Многие из них подтверждаются фактом, что экспериментальные данные не выявляют однозначную зависимость зонной структуры от конфигураций нанотрубок. Причиной тому может быть как технологическое несовершенство синтеза и анализа свойств 13

Рис. 3.3: Схема определения индексов хиральности материала, так и пренебрежение учётом межэлектронных взаимодействий в построении моделей. Высказывалось предположение, что к неверной интерпретации зависимости электропроводящих свойств нанотрубок от соотношения индексов хиральности могло привести пренебрежение в расчётах кулоновским взаимодействием внутри атома углерода [14]. Различают нанотрубки, обладающие зеркальной симметрией — типа «зигзаг» (англ. zigzag) с индексами (n, 0) (рис. 3.4; свёртка под углом 0◦ ) и типа «кресло» (англ. armchair) с индексами (n, n) (рис. 3.5; свёртка под углом 30◦ ) — их также именуют ахиральными или прямыми. Остальные наборы определяют не зеркально симметричные, то есть хиральные (спиральные) нанотрубки (рис. 3.6). Также нанотрубки подразделяют на открытые и закрытые. В последнем случае на торцах отсутствуют свободные связи. Ранее речь шла об однослойных нанотрубках. Первоначально же были получены многослойные (с двумя, пятью и семью концентрическими стенками) нанотрубки. Их структура разнообразна и зависит от условий синтеза. Обычно говорят о двух основных видах многослойных нанотру14

Рис. 3.4: Пример прямой нанотрубки типа «зигзаг» (5, 0) Рис. 3.5: Пример прямой нанотрубки типа «кресло» (3, 3) Рис. 3.6: Пример спиральной нанотрубки (5, 3) бок: о представляющих собой несколько вложенных одностенных или являющихся свёрнутым в свиток графеновым листом (рис. 3.7). Точной классификации структур многослойных нанотрубок не существует, так как изображения боковых сечений неоднозначны, а получение изображений торцевых срезов затруднительно. Известно несколько методов получения углеродных нанотрубок [5]. Изначально их получали электродуговым способом (С. Ииджима), подобно фуллеренам, что приводило к смесям однослойных с многослойными. Затем был предложен метод лазерной абляции графита в присутствии частиц металла, выступающих в качестве катализатора. Этот способ позволяет получать одностенные нанотрубки. В последнее время активно развиваются 15

Рис. 3.7: Формы многослойных нанотрубок подходы, основанные на осаждении из газовой фазы — они считаются наиболее коммерчески перспективными. Рассмотрим некоторые вопросы, касающиеся особенностей углеродных нанотрубок, используемых в качестве полевых эмиттеров. Многие исследователи сходятся во мнении, что поведение эмиттеров на основе углеродных нанотрубок во внешнем электрическом поле заметно отличается от поведения металлов. Одна из причин может быть связана с толщиной нанотрубок. Поскольку стенки материала очень тонкие, электрическое поле может проникать в трубку и изменять традиционный механизм автоэмиссии [15]. Кроме того, для реальных нанотрубок с несколькими тысячами электронов из-за ограниченности вычислительных ресурсов становится проблематичным моделирование целого объекта. Группой китайских учёных в таком случае было предложено рассматривать наконечник нанотрубки отдельно [16]. Целью этих исследований является нахождение способов получения наибольшего эмиссионного тока. На сегодняшний день известно, что к этому приводит добавление атомов некоторых элементов [17]. 3.2.2. Карбидокремниевые нанотрубки Из-за бесчисленных потенциальных приложений углеродные нанотрубки, несомненно, представляют большой научный и коммерческий интерес. Однако одним из ограничений их использования является разрушение материала в высокотемпературных условиях и агрессивных средах. Именно для таких случаев разрабатываются карбидокремниевые нанотрубки. К 16

примеру, последние остаются стабильными к окислению до 1000 ◦ С, в то время как первые — только до 600 ◦ C. Карбид кремния с атомистической точки зрения представляет связующее звено между кремниевой и углеродной электроникой. Из-за большой энергии связи кремния с углеродом этот материал сложен в изготовлении и обработке, но более долговечен. Нанотрубки на основе карбида кремния обладают теми же необычными свойствами и преимуществами углеродных нанотрубок перед другими эмиссионными материалами, что возникают вследствие атомной структуры, а именно из-за высокого отношения длины модели к диаметру. Визуально карбидокремниевые нанотрубки отличается от углеродных только тем, что места половины атомов в них занимают атомы кремния (рис. 3.8), однако их структура принципиально иная и имеет не грефеноподобную гибритизацию sp2 , а кристаллообразную — sp3 . В остальном для них применима та же терминология, включая индексы хиральности. Рис. 3.8: Закрытая нанотрубка карбида кремния (3, 3) Карбид кремния традиционно используется в качестве конструкционных материалов и компонентов керамики, но также представляет интерес для разработки новых компонентов высокотемпературных электронных устройств, так как обладает высокой термической и коррозийной устойчивостью в сочетании с большой величиной запрещённой зоны [18]. В последнее время резко возрос интерес к наноразмерным структурам карбида кремния. Одним из направлений здесь являются работы, посвящённые возможности создания плёнки графена на поверхности катода методом термического разложения [19], [20]. Это помогает уменьшить 17

негативное воздействие разрушающих факторов, увеличить долговечность катода и получаемого эмиссионного тока. Теоретические исследования в этой области проводятся в основном с целью определения равновесных структур. Это связано, в первую очередь, с разнообразием политипных разновидностей карбида кремния и стремлением к их классификации [2]. Как правило, выделяют три структурных семейства карбидокремниевых фаз: • листы карбида кремния, подобные графену (L-фазы, от англ. layer — слой), • однослойные нанотрубки из карбида кремния (T-фазы, от англ. tube — трубка), • фуллереноподобные молекулы (C-фазы, от англ. cluster — кластер). Особое место здесь занимают не нанотрубки, а фулсицены — шарообразные структуры карбида кремния, подобные фуллеренам [21]. Однако с точки зрения эмиссионной электроники наиболее перспективными являются именно нанотрубки. Существует несколько методик получения нанотрубок карбида кремния. Среди них можно выделить химическое превращение углеродных нанотрубок в карбидокремниевые, прямое выращивание на катализаторе и синтез из твердотельных пленок [2]. Для реализации последнего используют два основных подхода, суть которых состоит в освобождении тонкой пленки от подложки с помощью селективного травления. 18

4. Методы исследования 4.1. Математическая модель задачи В квантовой механике состояние системы описывается волновой функцией, квадрат модуля которой имеет смысл плотности вероятности нахождения микрочастицы в некоторой точке пространства. Волновая функция удовлетворяет уравнению Шрёдингера — дифференциальному уравнению в частных производных. Нерелятивистское уравнение Шрёдингера в координатном представлении для точечной частицы массы m, движущейся в поле с потенциалом U (r, t), в атомных единицах имеет вид: − ∂ψ(r, t) 1 ∆ψ(r, t) + U (r, t)ψ(r, t) = i , 2m ∂t (4.1) где ψ(r, t) — волновая функция системы, i — мнимая единица, ∆ — оператор Лапласа. Одну из главных характеристик автоэлектронной эмиссии — плотность эмиссионного тока — можно найти через волновую функцию. Покажем это. Как известно, волновая функция изменяется в пространстве и времени в соответствии с законом сохранения [22]. Сформулируем его, найдя производную по времени от плотности вероятности. Для этого запишем комплексно сопряжённое к (4.1) уравнение 1 ∂ψ ∗ ∗ ∗ − ∆ψ + U (r, t)ψ = −i , 2m ∂t (4.2) и вычтем из (4.1) (4.2), предварительно умножив на ψ ∗ и ψ соответственно: − ∂ψ 1 ∗ ψ ∆ψ + U (r, t)ψψ ∗ = iψ ∗ , 2m ∂t 1 ∂ψ ∗ ∗ ∗ − ψ∆ψ + U (r, t)ψψ = −iψ . 2m ∂t Получим ∗ 1 ∂ψ ∂ψ − (ψ ∗ ∆ψ − ψ∆ψ ∗ ) = i ψ − ψ∗ , 2m ∂t ∂t 19

или, переписывая, 1 ∂ψ ∗ ψ ∗ ∗ ∇ (ψ ∇ψ − ψ∇ψ ) + = 0, 2mi ∂t откуда, обозначив выражение в скобках за j= 1 (ψ ∗ ∇ψ − ψ∇ψ ∗ ), 2mi выводим уравнение непрерывности: ∂|ψ|2 div j + = 0. ∂t Последнее выражение имеет форму закона сохранения плотности вероятности, где div j — удельная мощность источников в поле вектора j. Вектор j есть плотность потока вероятности, или, по аналогии с законом сохранения, плотность электрического тока, поскольку масса и заряд частицы неразделимы. Отсюда имеем, что связь эмиссионного тока с волновой функцией может быть записана в виде j= 1 (ψ ∗ ∇ψ − ψ∇ψ ∗ ). 2mi (4.3) Таким образом, математически наша задача состоит в интегрировании уравнения Шрёдингера с целью нахождения волновой функции. Обозначим H=− 1 ∆ + U (r, t) 2m и перепишем (4.1) в виде Hψ(r, t) = i ∂ψ(r, t) . ∂t (4.4) Оператор полной энергии системы H называется оператором Гамильтона (гамильтонианом). Предположим, что гамильтониан системы не зависит от времени, и разделим переменные ψ(r, t) = Ψ(r)Φ(t). Тогда (4.4) примет вид Φ(t)HΨ(r) = iΨ(r) 20 ∂Φ(t) , ∂t

или, иначе 1 ∂Φ(t) 1 HΨ(r) = i . Ψ(r) Φ(t) ∂t (4.5) Левая часть (4.5) является функцией только координат, а правая – только времени. Поэтому обе части должны быть равны некоторой (одной и той же) постоянной, которую принято обозначать E. Разрешение уравнения относительно t даёт Φ(t) = e−iEt , следовательно, ψ(r, t) = e−iEt Ψ(r), а не зависящее от времени уравнение Шрёдингера имеет вид HΨ(r) = EΨ(r), (4.6) причём E имеет смысл полной энергии. Здесь и далее мы, не умаляя общности, будем рассматривать только стационарное уравнение Шрёдингера (4.6). Для системы из n электронов и p ядер оператор Гамильтона в атомных единицах имеет вид n p n p n n p p X X zj X X 1 X X zj zl 1X 1X H=− ∆i − ∆j − + + . (4.7) 2 i=1 2 j=1 r r r ij ik jl i=1 j=1 i=1 j=1 k>i l>j Первые две суммы отвечают кинетической энергии электронов и ядер, остальные три — потенциальным энергиям взаимодействия электронов с ядрами, электронов с электронами и ядер с ядрами соответственно (здесь z — зарядовое число). Сложный вид (4.7) затрудняет решение уравнения (4.6) уже для систем с более чем одни электроном. Нанотрубки, как и многие другие реальные системы, являются многоэлектронными. В таких случаях уравнение Шрёдингера решается с использованием приближённых методов. 21

4.2. Обзор методов решения уравнения Шрёдингера Можно предложить следующую классификацию приближённых методов решения уравнения Шрёдингера (таб. 4.1). Таблица 4.1: Классификация приближённых методов решения уравнения Шрёдингера По выбору математического аппарата Методы теории возмущений Вариационные методы По привлечению экспериментальных данных Неэмпирические Полуэмпирические По способу формирования волновой функции Метод валентных связей Метод молекулярных орбиталей Остановимся на каждом из них. Методы теории возмущений состоят в нахождении решения в виде ряда по степеням члена, именуемого возмущением, где главный член ряда соответствует точно решаемой задаче, а остаточный описывает погрешность приближения. Этот метод эффективен в случае сложного гамильтониана, близкого к тому, что может быть найден точно, и неприменим, если имеется особенность в точке разложения. Метод функционала электронной плотности, используемый в работе, базируется на ином подходе — вариационном — и состоит в нахождении решения путём минимизации функционала (в данном случае, функционала энергии). Вариационным методам в квантовой механике будет посвящён один из последующих разделов. Как будет показано позднее, гамильтониан системы можно заметно упростить исходя из естественных предположений, используя приближение Борна—Оппенгеймера. Такой подход характерен для всех квантовохимических методов, однако некоторые из них этим упрощением не ограничиваются. Отсюда возникает разделение методов на неэмпирические (не использующие иные, кроме оговоренного, приближения) и полуэмпирические. Последние включают ряд других приближений, связанных в основ22

ном с валентным приближением (описанием движения только валентных электронов) или пренебрежением определенными классами молекулярных интегралов в точных уравнениях системы, а также с использованием данных, полученных опытным путём. Полуэмпирические методы позволяют ускорить расчёты ценой точности решения. Обзор этих методов выходит за рамки данной работы и подробно рассмотрен не будет. Мы остановимся только на неэмпирических методах исследования структур, называемых также методами ab initio или, иначе, методами из первых принципов. Различие методов приближённого решения уравнения Шрёдингера связано также с формами химического моделирования. В первой половине XX века объединение результатов частных решений уравнения Шрёдингера, экспериментальных данных, плодов изучения молекулы водорода и распространения на случаи более сложных молекул позволили установить, что образование химической связи является результатом взаимопроникновения электронных облаков — процессом спаривания электронов с противоположно направленными спинами [23]. Такое представление и схематичное изображение легло в основу метода валентных связей. Обладая наглядностью, информацией о способности к образованию ковалентных связей, её направленности и структуре молекулы, в ряде случаев он приводит к неверным заключениям. Другой, более общий способ описания – метод молекулярных орбиталей. Полагается, что состояние электронов в атоме представляет собой совокупность атомных электронных облаков (орбиталей), а каждая орбиталь характеризуется набором четырёх квантовых чисел (главным, азимутальным, магнитным и спиновым). Метод молекулярных орбиталей исходит из предположения, что состояние электронов в молекуле (аналогично атомному) может быть описано совокупностью молекулярных электронных орбиталей со своим набором молекулярных квантовых чисел. Молекулярное облако может быть децентрализовано: сосредоточено вблизи одного из атомных ядер; может принадлежать нескольким атомным ядрам и образовывать многоцентровую химическую связь. Это объясняет широту применения метода молекулярных орбиталей в сравнении с методом валентных связей, опирающемся на двуцентровую связь как частных 23

случай многоцентровой. 4.3. Математический аппарат 4.3.1. Вариационный принцип квантовой механики Задачу решения уравнения Шрёдингера, являющегося дифференциальным, можно свести к задаче на экстремум функционала, которым и является среднее значение оператора Гамильтона (энергии) на классе функций, имеющих смысл волновых (непрерывных с первыми производными, ограниченных и интегрируемых с квадратом). Вариационный принцип квантовой механики можно сформулировать следующим образом: если наименьшее значение гамильтониана системы равно E0 , а ψ0 — соответствующая собственная функция (точная волновая функция), то для любой произвольной функции Φ, называемой пробной, выполняется условие: Z E= Z ∗ Φ HΦdV ≥ ψ0∗ Hψ0 dV = E0 , причем точное равенство E = E0 достигается, если пробная функция совпадает с точной волновой Φ = ψ. Действительно, пусть функция Φ представима в виде линейной комбинации функций χi : Φ= N X c i χi , i=1 где χi удовлетворяют уравнению Hχi = Eχi . Найдём среднее значение энергии для пробной функции Φ: Z E= ∗ Φ HΦdV = X Z cn ck χ∗n Hχk dV = n,k = X n,k 24 Z cn ck Ek χ∗n χk dV = X n |cn |2 En .

Запишем условие нормировки пробной функции: Z Z X X X ∗ |cn |2 = 1, cn ck χ∗n χk dV = cn ck δnk = Φ ΦdV = n,k n n,k то есть X |cn |2 = 1. n Имеем E= X |cn |2 En , E0 = X |cn |2 E0 , n n вычитая, получаем E − E0 = X (En − E0 )|cn |2 . n Так как E0 соответствует основному состоянию и, следовательно, имеет минимальное значение среди всех энергий En , то верно En ≥ E0 . Поскольку |cn |2 ≥ 0 всегда, то и E − E0 ≥ 0. Именно использование пробной функции позволяет решить уравнение Шрёдингера, если это невозможно проделать аналитически; точность решения определяется точностью выбора пробной функции. В общем случае пробная функция зависит от вариационных коэффициентов. 4.3.2. Вариационный метод Ритца Существует два основных пути решения вариационных задач: сведение к решению дифференциальных уравнений или использование прямых методов по поиску функции, доставляющей экстремальное значение функционалу на заданном классе. К последним относится наиболее популярный метод Ритца, состоящий в построении минимизирующей последовательности в виде суперпозиции некоторых известных функций, удовлетворяющих тем же условиям, что и предполагаемое решение [24]. Будем искать пробную волновую функцию в виде линейной комбинации (линейный вариационный метод): ψi = N X j=1 25 Cij χi ,

где χi — некоторые известные независимые функции (базисные). Такое представление оправдано тем, что оператор Гамильтона является линейным самосопряжённым и имеет полную систему собственных функций P∞ Ψ = n=1 Cn ψn . Необходимо найти такие коэффициенты Cij , которым соответствует наименьшее значение энергии. Если функция Φ не нормирована, а χn и χk не ортогональны, то среднее значение энергии равно R ∗ P P R ∗ c c Hχ dV χ Φ HΦdV n k k n,k cn ck Hnk n,k R n P = P = , E= R ∗ Φ ΦdV χ∗n χk dV n,k cn ck Snk n,k cn ck R где Hnk = χ∗n Hχk dV — матричные элементы гамильтониана, Snk = R ∗ χn χk dV — матричные элементы интегралов перекрывания функций. Переписывая это выражение в виде E X cn ck Snk = n,k X cn ck Hnk n,k и дифференцируя по коэффициентам cn , получаем X ∂E n,k ∂cn cn ck Snk + E X ck Snk = k X ck Hnk . k Так как производная энергии по вариационному коэффициенту в точке экстремума равна нулю, окончательно записываем X ck (Hnk − ESnk ) = 0, n, k = 1, . . . , N. (4.8) k Чтобы это уравнение имело нетривиальное решение, необходимо и достаточно выполнения условия |Hnk − ESnk | = 0. Наименьший из корней этого уравнения, называемого секулярным или вековым, соответствует энергии основного состояния, остальные корни представляют энергии возбуждённых состояний. Если пробная функция подобрана удачно и варьируемые параметры обеспечивают достаточную гибкость функции Ψ, то улучшение пробной функции может привести к точному решению уравнения Шрёдингера. 26

4.4. Особенности теоретических методов исследования Выше мы уже отметили преимущества и недостатки неэмпирических методов перед полуэмпирическими и частично обосновали необходимость использования первых. Несмотря на трудоёмкость, методы из первых принципов являются более последовательными и общими, а также обеспечивающими более точные результаты [11]. Этот раздел предваряет описание теоретических основ данного класса методов и призван выявить их применимость в задачах моделирования поведения материала эмиттера в условиях сильных электрических полей. Разумеется, предсказание свойств веществ на основе знания только их химического состава давно привлекает исследователей. Теоретический фундамент для этого был создан столетие назад, но только значительные успехи, достигнутые в развитии вычислительной техники и методах решения больших систем уравнений, позволили в настоящее время получить при расчёте свойств веществ из первых принципов точность, сопоставимую с экспериментальной. Ab initio методы позволяют прогнозировать поведение веществ в экстремальных условиях (например, при сверхвысоких давлениях, которые пока недостижимы в лабораторях), изучать свойства опасных (радиоактивных, взрывоопасных) веществ. Они чрезвычайно полезны для прогнозирования свойств новых соединений, которые пока ещё не синтезированы. Кроме того, расчёты из первых принципов позволяют глубже разобраться в сути физических явлений, происходящих в уже исследованных материалах. Главным достоинством расчётов из первых принципов является описание атомного взаимодействия с учётом квантовых эффектов. Несмотря на значительный прогресс в моделировании из первых принципов со времён ранних работ, датируемыми последними двадцатью годами, количественное описание процессов с их помощью и по сей день остаётся проблемой. Известно, что моделирование методом функционала электронной плотности, на котором основаны настоящие исследования, неплохо описывает качественную картину, но количественно оказывается 27

далёким от совершенства [8]. С использованием методов ab initio для наноэмиттеров связаны две основные сложности. Во-первых, невозможность пренебрежения атомной структурой эмиттера по сравнению с макроскопической. Это легко проиллюстрировать следующим примером: в классической модели полевой электронной эмиссии коэффициент усиления поля и работа выхода полагаются постоянными. На самом деле они зависят от напряжённости поля. Кроме того, поскольку имеют место множественные межэлектронные взаимодействия, традиционная картина одного туннелирующего сквозь потенциальный барьер электрона не приемлема. Более точное описание процесса обеспечивает концепция квазичастиц, которая может быть удобно представлена на языке теории функционала электронной плотности. Приступим к обзору основных вычислительных методов квантовой химии. 4.5. Вычислительные методы квантовой химии 4.5.1. Метод Хартри Как отмечалось ранее, сложный вид оператора Гамильтона многоэлектронной системы является препятствием к решению уравнения Шрёдингера и вынуждает использовать упрощения. Наиболее общим для анализа молекулярных систем является приближение Борна—Оппенгеймера, которое состоит в том, что описание ядер атомов и электронов в системе происходит отдельно. Это естественное предположение, вытекающее из значительной разницы масс, а, следовательно, и скорости движения. Таким образом, систему можно рассматривать с точки зрения движения электронов в электростатическом поле покоящихся ядер. В этом случае оператор Гамильтона для системы из n электронов и p ядер n p n p n n p p X X zj X X 1 X X zj zl 1X 1X ∆i − ∆j − + + H=− 2 i=1 2 j=1 r rik j=1 rjl i=1 j=1 ij i=1 k>i 28 l>j

упростится следующим образом: n n p n n X X zj X X 1 1X Hel = − ∆i − + , 2 i=1 r r ij ik i=1 j=1 i=1 k>i Hnuc = p X p X zj zl j=1 l>j rjl . Идея метода самосогласованного поля, предложенная английским учёным Д. Хартри, состоит в замене межэлектронного взаимодействия моделью, согласно которой каждый электрон движется в некотором усреднённом поле, создаваемом другими электронами. Это позволяет упростить ещё одно слагаемое в операторе Гамильтона — потенциал межэлектронноP P го взаимодействия ni=1 nk>i r1ik — и заменить его выражением, зависящим от координат только одного электрона. Волновую функцию системы (полную) будем искать в виде произведения волновых функций, соответствующих каждому электрону: ψ(r1 , . . . , rn ) = ψ(r1 ) . . . ψ(rn ). Такая форма записи предполагает, что электрон в атоме движется независимо от других, но поскольку это не так, погрешность расчёта полной энергии, обусловленная пренебрежением электронной корреляцией, должна учитываться в дополнительном слагаемом в гамильтониане. Более того, такое представление нарушает принцип запрета Паули. Последний недостаток устраняется в методе Хартри—Фока. 4.5.2. Метод Хартри—Фока Идея метода аналогична методу Хартри с той разницей, что в нём учтена антисимметрия волновой функции. Это возможно благодаря представлению её в виде определителя Слэтера:

ψ1 (r1 ) ψ1 (rn )

ψ(r1 , . . . , rn ) =

. . . . . . . . . . . . . . .

ψn (r1 ) ψn (rn )

Действуя 2n-электронным гамильтонианом на волновую функцию в виде детерминанта Слэтера, получим выражение для полной энергии Z E= ∗ ψ HψdV = n X Hi + n X n X i=1 (2Jik − Kij ), j=1 i>j где Hi — остовный интеграл, являющийся суммой кинетической энергии электрона на орбитали ψi и потенциальной энергии притяжения его к ядру, удвоенный, так как на каждой орбитали располагается два электрона: n z 1X ∆i − ; Hi = − 2 i=1 ri Jij — кулоновский интеграл, представляющий среднюю энергию электростатического отталкивания электронов, находящихся на орбиталях ψi и ψj (удвоенный по той же причине): Jij = 1 ; ri j Kij — обменный интеграл, частично учитывающий электронную корреляцию: он определяет величину, на которую в меньшую сторону изменяется электростатическая энергия отталкивания электронов вследствие принципа запрета Паули. Но по этой же причине и среднее расстояние между электронами больше на величину обменного интеграла. В методе Хартри— Фока корреляция, вызванная кулоновским отталкиванием пар электронов с противоположно направленными спинами, не учитывается. Орбитали определяются из условия минимума полной энергии вариационным методом. Составим вспомогательную функцию F =E− n X n X Z εij ψi ψj dV = E − i=1 j=1 n X Z εi ψi2 dV . i=1 Условие минимума — обращение в нуль первой вариации: δF = 2 n X i=1 δHi + n X n X (2δJij − δKij ) − j=1 i>j n X i=1 30 Z εi δ(ψi2 )dV = 0.

С использованием выражений для кулоновского и обменного интегралов получаем n Z X  δψi 2Hi (1)ψi (1) + i=1 n X Z 2ψi (1) j6=i ψj2 (2) r12 Z dV2 − ψj (1) ψi (2)ψj (1) dV2 r12 где символы 1 и 2 обозначают набор аргументов. Для любых ψi это выражение справедливо, если в нуль обращается выражение в фигурных скобках: 2Hi (1)ψi (1) + n X Z 2ψi (1) ψj2 (2) r12 j6=i ! Z ψi (2)ψj (1) dV2 − r12 dV2 − ψj (1) − εi ψi (1) = 0, i = 1, . . . , n. Эта система уравнений называется Хартри—Фоковской и решается итерационным путём. Каждое уравнение содержит координаты только одного электрона, но нуждается в указании эффективного потенциала межэлектронного отталкивания 2Jij − Kij , зависящего от функций ψj . Схему решения можно представить следующим образом: ψj0 Jj0 , Kj0 , 2Jj0 − Kj0 Jjk , Kjk , ψjk 2Jjk − Kjk , k := k + 1 ψjk ≈ ψjk−1 Нет Да ψjk 31 ! − εi ψ

Основной недостаток метода Хартри—Фока состоит в том, что функции ψi не имеют аналитического вида и могут быть получены только в виде дискретных значений, а обменная энергия учитывается частично. Значения εi = Hi + n X (2δJij − δKij ) j=1 соответствуют полной энергии электрона, находящегося на орбитали ψi . Выражение для полной энергии системы имеет вид E=2 n X i=1 εi − n X n X (2Jij − Kij ). j=1 i>j 4.5.3. Метод функционала электронной плотности Резюмируем предыдущие результаты. Метод самосогласованного поля приводит к замене одноэлектронного уравнения Шрёдингера на уравнение движения электрона в самосогласованном потенциале: ∆ − + Vef f (r) ψi = εi ψi . 2 В приближении Хартри, игнорирующем корреляционные эффекты и фермионную природу электронов (отвечающие за дополнительное притяжение), этот потенциал определялся как сумма электростатических взаимодействий данного электрона с другими электронами (отталкивание) и со всеми ядрами окружающих атомов (притяжение): Z ρ(r0 ) Vef f (r) = Vex (r) + dr0 . 0 |r − r | В приближении Хартри—Фока обменный член явно включён в уравнения. Это частично решает проблему некорректного физического моделирования, но изменяет вид уравнений так, что они перестают быть линейными относительно ψi . Вследствие этого возникает потребность в более мощных, чем для метода Хартри, алгоритмах диагонализации систем линейных уравнений, в итоге — к высоким вычислительным затратам. Потому метод Хартри—Фока целесообразно применять только для небольших 32

молекул (вычислительная сложность метода составляет N 4 , где N — число базисных функций или атомов). Необходимость в создании подхода, учитывающего все указанные физические явления и не приводящего к нелинейным уравнениям, привела к созданию метода функционала электронной плотности. Первой теорией, показавшей путь решения данной проблемы, явилась теория Томаса— Ферми, не учитывающая обмен, но позволяющая найти кинетический вклад в гамильтониан без явного решения уравнения Шрёдингера. Теория Томаса—Ферми Допустим, что моделью любой молекулярной или кристаллической структуры можно считать электронную жидкость с включенными в неё ядрами атомов. Пусть электроны полностью делокализованы, то есть их волновые функции являются комбинацией плоских волн, и справедлив подход самосогласованного поля (уравнения движения электрона во внешнем потенциале справедливы для потенциалов ядер в виде δ-функций). Плотность однородного электронного газа имеет вид ρ= 1 (2εf )3 /2, 2 3π плотность кинетической энергии системы невзаимодействующих электронов равна 1 t[ρ(r)] = (2π)3 kf Z 0 3 k 2 dk = ρεf . 5 Таким образом, кинетическая энергия может быть представлена функционалом Z t[ρ(r)] = Z t[ρ(r)]dr = C ρ5/3 (r)dr, C= 3 (3π 2 )2/3 . 10 Функционал энергии, связанный с межэлектронным взаимодействием (вклад электростатического кулоновского взаимодействия Хартри), имеет вид 1 Eee = 2 Z Z ρ(r)ρ(r0 ) drdr0 . 0 |r − r | В предположении о движении электронов во внешнем потенциале его 33

энергия представима в виде Z Eex = ρ(r)Vex (r)dr. Минимизация функционала полной энергии системы E[ρ] при условии соR хранения числа электронов N = ρ(r)dr методом неопределённых множителей Лагранжа Z E[ρ(r)] = T [ρ(r)] + Eee [ρ(r)] + Z ρ(r)Vex (r)dr + λ N − ρ(r)dr даёт уравнение Томаса—Ферми Z ρ(r0 ) 5 2/3 Cρ (r) + dr0 + Vex (r) − λ = 0, 3 |r − r0 | (4.9) где λ ассоциируют с химическим потенциалом. 3 3/2 ) и Vef f (r) = Vex (r) + Через обозначения γ = ( 5C R ρ(r0 ) 0 |r−r0 | dr элек- тронная плотность выражается в виде ρ(r) = γ(λ − Vef f (r))3/2 . Таким образом было установлено, что между электронной плотностью и внешним потенциалом существует взаимно однозначное соответствие. Затем этот факт был теоретически доказан учёными В. Коном и П. Хоэнбергом. Позднее также было показано, что теория Томаса—Ферми является точной в случае предела бесконечных зарядов ядер [25]. Для реальных систем она полезна в описании некоторых качественных характеристик (например, полных энергий атомов), но её количественные предсказания не воспроизводят химические и материаловедческие особенности строения: теория неприменима к валентным электронам и не свидетельствует о наличии химической связи. Важным в таком подходе является то, что выражение (4.9) сформулировано в терминах электронной плотности, а не волновых функций. Этот факт стал предтечей метода функционала электронной плотности. 34

Итак, электронная плотность может быть получена через волновую функцию интегрированием по всем пространственным координатам, кроме одной: Z ρ(r) = n Z ... |ψ(r1 , . . . , rn )|2 dr2 . . . drn . Эта величина в отличие от волновой функции имеет физический смысл (то есть может быть измерена) и зависит только от трёх пространственных переменных для любой n-частичной системы, у которой волновая функция в общем случае зависит от 4n координат. Следующим важным шагом в развитии теории функционала электронной плотности явилось установление двух теорем Кона—Хоэнберга в 1964 г. [26]. Сформулируем их. Первая теорема Кона—Хоэнберга Для любой системы взаимодействующих электронов, находящихся во внешнем потенциале, этот потенциал определяется однозначно (с точностью до произвольной постоянной) электронной плотностью основного состояния. Доказательство. I Пусть ρ(r) — плотность основного невырожденного состояния ψ1 с энергией E1 системы n электронов в потенциале Vex1 (r). Тогда по определению энергии Z Z Z E1 = ψ1∗ H1 ψ1 dV = Vex1 (r)ρ(r)dr + ψ1∗ (T + Π)ψ1 dV, где n n X n X 1 Π= . r ij j=1 i>j 1X ∆i , T =− 2 i=1 Пусть существует второй потенциал Vex2 (r) 6= Vex1 (r) + c, c = const, отвечающий состоянию ψ2 6= ψ1 с плотностью ρ(r). Тогда Z Z Z E2 = ψ2∗ H2 ψ2 dV = Vex2 (r)dr + ψ2∗ (T + Π)ψ2 dV. 35

Если состояние ψ1 является невырожденным, по вариационному принципу получим: Z Z Z E1 < ψ2∗ H1 ψ2 dV = Vex1 (r)dr + ψ2∗ (T + Π)ψ2 dV = Z = E2 + (Vex1 (r) − Vex2 (r))ρ(r)dr. Аналогично Z Z Z E2 ≤ ψ1∗ H2 ψ1 dV = Vex1 (r)dr + ψ2∗ (T + Π)ψ2 dV = Z = E1 + (Vex1 (r) − Vex2 (r))ρ(r)dr (знак нестрогий, так как невырожденность состояния ψ2 не предполагается). Сложение этих неравенств приводит к противоречию E1 +E2 < E1 +E2 , что означает ошибочность предположения и доказывает взаимную однозначность соответствия электронной плотности волновой функции. J Замечание: в случае вырожденного состояния теорема верна для плотности ρ(r) любого основного состояния. Вторая теорема Кона—Хоэнберга Существует универсальный функционал электронной плотности, справедливый для любого внешнего потенциала, причём в случае определённого потенциала экстремум достигается для электронной плотности основного состояния. Теорема по сути представляет собой вариационный принцип квантовой механики, сформулированный для функционала плотности, поэтому приводится без доказательства. Замечание: теоремы сформулированы для основного состояния электронной подсистемы в отсутствие магнитного поля, но могут быть обобщены путём введения зависимости от времени — это позволяет рассчитывать состояния возбужденных электронов [27], [28]. Таким образом, функционал полной энергии может быть записан в виде Z E[ρ(r)] = T [ρ(r)] + Eee [ρ(r)] + 36 ρ(r)Vex (r)dr.

Часть слагаемых T [ρ(r)] + Eee [ρ(r)] = FHK [ρ(r)] именуют функционалом Кона—Хоэнберга, вид которого не зависит ни от числа частиц, ни от внешнего потенциала. Знание его точного значения позволяет получить точную электронную плотность основного состояния, но на практике достичь этого невозможно. Теоремы гласят, что энергия основного состояния является функционалом электронной плотности и достигает минимума в случае её точного значения, но не дают способа построения этого функционала. Движение в направлении практической реализации метода было предложено учёными В. Коном и Л. Шэмом. Уравнения Кона—Шэма В 1965 г. В. Кон и Л. Шэм показали, что плотность основного состояния системы взаимодействующих частиц может быть вычислена как плотность в основном состоянии вспомогательной невзаимодействующей системы [29]. Иначе говоря, для каждой взаимодействующей системы с внешним потенциалом существует локальный потенциал, такой, что плотность системы невзаимодействующих частиц равна плотности взаимодействующей системы. Для вычисления этого потенциала необходимо знать функцию электронной плотности, которую можно получить из орбиталей. Но они, в свою очередь, ищутся через решение уравнения Шрёдингера, в котором гамильтониан завязан с потенциалом. Следовательно, требуется самосогласованное итерационное решение. В случае невзаимодействующих электронов FHK [ρ(r)] = T [ρ(r)] и энергию и плотность основного состояния можно найти, вычислив собственные функции и собственные значения εi , i = 1, . . . , n одночастич- ных уравнений вида ∆ − + Vex (r) ψi = εi ψi . 2 37

В результате E= n X εi , ρ(r) = i=1 n X |ψi (r)|2 i=1 и энергия основного состояния представима как n X Z εi = T [ρ(r)] + ρ(r)Vex (r)dr. i=1 В случае взаимодействующих электронов функционал Кона—Хоэнберга имеет вид FHK [ρ(r)] = T [ρ(r)] + Eee [ρ(r)] + Eex [ρ(r)], а система уравнений — ∆ − + Vef f (r) ψi = εi ψi . 2 Алгоритм метода функционала электронной плотности состоит из следующих шагов: • задание входных данных: координат атомов, зарядов ядер, полного числа электронов, формы псевдопотенциала (если используется); • задание начальной электронной плотности (как суммы зарядов заданной первоначальной атомной конфигурации или из предварительных полуэмпирических расчетов); • построение потенциала Хартри, например, с помощью решения уравнения Пуассона ∆VH (r) = −4πρ(r); • построение обменного и эффективного потенциалов; • решение уравнений Кона—Шэма; • проверка сходимости ρx−1 ≈ ρx (интегральная разница между полученной плотностью и плотностью на предыдущей итерации по модулю меньше некоторого изначально заданного значения), • повторение процедуры с улучшенной волновой функцией до достижения критерия сходимости. 38

Общую схему решения методом функционала электронной плотности можно представить следующим образом: Задание базисного набора χ1 , . . . χN Вычисление Hnk , Snk Формирование пробной функции PN ψn = j=1 Cnj χn Решение уравнений Кона—Шэма P k ck (Hnk − εn Snk ) = 0 n := n + 1 Нет Расчёт εn , E Сходимость n := n + 1; по энергии k := k + 1 Да Сходимость по координатам Да ρ 39 Нет

Приближения для обменно-корреляционного члена Обменно-корреляционный член в потенциале содержит всё, что не вошло в предыдущие, и теория верна тогда, когда он точно известен. Различные формы методов семейства функционала электронной плотности различны только в его выборе [30]. На практике приходится полагаться на некоторое приближённое значение обменно-корреляционной энергии. Как правило, её ищут в виде Z Exc [ρ] = Ex [ρ] + Ec [ρ] = Z ρ(r)εx (ρ(r))dr + ρ(r)εc (ρ(r))dr. Перечислим некоторые известные приближения для обменно-корреляционного функционала. • Приближение локальной плотности (LDA) Z LDA Exc [ρ(r)] = ε[ρ(r)]ρ(r)dr. Здесь ε[ρ(r)] представляет обменно-корреляционную энергия на единицу объема гомогенного электронного газа плотности ρ(r). Предположение об однородности является приближением, поскольку внутри молекулы заряды распределены неравномерно. • Приближение локальной спиновой плотности (LSDA) — обобщение предыдущего LDA приближения εLSDA [ρ(r)] = −21/3 Cx ρ1/3 x α + 1/3 ρβ , 1/3 3 3 − . Cx = 4 π • Обобщенное градиентное приближение (GGA) связывает Ex и Ec не только с плотностью ρ(r), но и с её первой и второй производными. – Обменная поправка Бекке (Becke, 1988): LDA εB + ∆εB x = εx x, x2 |∇ρ| = −βρ , x = , 1 + 6βxsh−1 x ρ4/3 параметр β подбирается по известным данным для атомов. ∆εB x 1/3 40

– Поправка Пердью и Ванга (1991) обменная: 2 εPx W = εLDA x 1 + xa1 sh−1 (xa2 ) + (a3 + a4 e−bx x2 ) 1 + xa1 sh−1 (xa2 ) + a5 x2 где a1 , a2 , a3 , a4 , a5 , b — подобранные константы, x = ! , |∇ρ| ; ρ4/3 корреляционная: εP W = εLDA + ∆εP W . – Градиентная коррекция Ли, Янга и Парра: ρα ρ α P − abω∗ εLY = −4a 2 c ρ (1 + dρ−1/3 ) ρ α ρα 2 2 ∗ X + ρ |∇ρα |2 + |∇ρβ |2 − |∇ρ|2 − ρ2α |∇ρβ |2 + (ρ2β |∇ρα |2 , 18 3 где 8/3 2 X = 144(22/3 )CF (ρ8/3 α + ρβ ) + (47 − 7δ)|∇ρ| − − (45 − δ) |∇ρα |2 + |∇ρβ |2 + 2ρ−1 ρα |∇ρα |2 + ρβ |∇ρβ |2 , −1/3 e−cρ dρ−1/3 3 1/3 (3π 2 )2/3 , ω= , δ = cρ + , C = F −1/3 14/3 −1/3 10 (1 + dρ )ρ 1 + dρ a, b, c, d — постоянные параметры, полученные экспериментально для атома гелия. • Гибридные методы hyb LDA Exc = Exc +a0 (ExHF −ExLDA )+ax (ExGGA −ExLDA )+ac (EcGGA −EcLDA ). Параметры a0 , ax , ac подбираются из условия лучшего воспроизведения термохимических свойств набора молекул. Наибольшее распространение среди гибридных методов получил трехпараметрический функционал Бекке: B3 Exc = (1 − a)ExLSDA + aExHF + b∆Ex B + EcV W N + c∆EcGGA . Последний член суммы представляет градиентную коррекцию корреляционного функционала. В зависимости от вида поправки получается 41

конкретная реализация гибридного метода. Например, использование поправки Пердью и Ванга даёт метод B3PW91; корреляционный функционал Ли, Янга и Парра приводит к B3LYP. Последний использован в данной работе и имеет вид B3LY P Exc = (1 − a)ExLSDA + aExHF + b∆Ex B + EcV W N + c∆EcLY P , где a, b, c — константы, подобранные Бекке экспериментально для набора простых химических соединений (G1-набор) [31]. Широкое практическое использование гибридных методов функционала электронной плотности обусловлено уникальным сочетанием низкой затратности приближения и высокой степени точности результатов расчётов. Альтернативные методы моделирования Как метод, являющийся расширением теории функционала электронной плотности на случай зависимых от времени систем, нельзя не отметить нестационарную теорию функционала электронной плотности. Подобно теоремам Кона—Хоэнберга, сформулированным для стационарного случая, существует теорема Рунге—Гросса, устанавливающая взаимно однозначное соответствие между зависящими от времени внешним потенциалом и электронной плотностью. Этот факт позволяет использовать аналогичный представленному аппарат. Поскольку для стационарного процесса при разделении переменных в уравнении Шрёдингера волновая функция имеет экспоненциальную зависимость от времени, учёт зависимости функционала энергии от времени в этом случае можно произвести разложением экcпоненты в ряд Тейлора. Иной подход к моделированию может обеспечить неравновесная теория функций Грина, менее популярная среди исследователей. Описание указанных методов выходит за рамки данной работы. 42

5. Практическая реализация 5.1. Формирование волновой функции Для начала вычислений (итерационного процесса) необходимо задать начальное приближение — пробную волновую функцию. Волновая функция системы или молекулярная орбиталь, говоря в терминах одноимённого метода, может быть представлена как сумма атомных орбиталей. Это оправдано практически, так как позволяет применить вариационный принцип Ритца, и логически следует из самого определения молекулярной орбитали. Такая форма метода молекулярных орбиталей называется методом линейной комбинации атомных орбиталей. Итак, пусть волновая функция системы представима в виде линейной комбинации волновых функций составляющих систему атомов — атомных орбиталей, представимых, в свою очередь, линейной комбинацией конечного числа базисных состояний χi : ψi = N X Cij χi , i = 1, . . . , n. (5.1) j=1 Именно выбор базисных атомных функций определяет точность такой аппроксимации. Поэтому базисные функции должны удовлетворять следующим требованиям: • наиболее точно описывать истинную волновую функцию и давать хорошее приближение возле ядер и на больших расстояниях от них, • допускать аналитическое вычисление интегралов, • не превышать некоторого количества, значительно затрудняющих вычисления или делающих их невозможными. Как правило, используют следующие наборы базисных функций: • плоские волны (∼ eikr ): – ортогонализованные плоские волны, – присоединённые плоские волны; 43

• обитали Слэтера (∼ e−βr ); 2 • орбитали Гаусса (∼ e−βr ); • численные орбитали, форма которых оптимизируется из атомных расчетов. 5.1.1. Плоские волны Только плоские волны формируют полный базисный набор, то есть с увеличением количества повышают точность решения. Разложение по плоским волнам удачно описывает электронную структуру периодичной кристаллической решётки, аморфных тел, атомных кластеров и систем с дефектами [28]. Выбор базиса в виде плоских волн является наиболее простым и естественным, но не обеспечивает быстрой сходимости ряда. Это связано с сильными осцилляциями волновых функций электронов проводимости в области центра иона, где они ведут себя как атомные волновые функции у валентных электронов. Для воспроизведения такого поведения необходимо использовать огромное число плоских волн (≈ 106 ). Но выбрав волновые функции, которые требуется найти, ортогональными к известным волновым функциям внутренних оболочек и обеспечив, тем самым, учёт осцилляции в центрах ионов, можно улучшить сходимость ряда [32]. Базис ортогонализованных плоских волн имеет вид:

1 1 opw − |φt (j)i φt (j)| √ , χk =

√ ikr Ωeikr Ωe t,j где φt (j) — нормированная собственная функция внутренних оболочек, Ω — нормированный объем кристалла, индекс t характеризует состояние внутренних оболочек (1s, 2s, 2p, . . . ), индекс j указывает на положение иона. Они ортогональны функциям внутренних оболочек в предположении, что волновые функции различных ионов не перекрываются (что является хорошей аппроксимацией), и образуют полную систему для разложения функций зоны проводимости. Присоединённые плоские волны в пространстве между сферами вы44

бираются как плоские волны, а в остальной области — как решения для сферически симметричного потенциала. Поэтому первоначально необходимо аппроксимировать потенциал, предположить его сферическую симметричность вблизи каждого ядра и малость радиуса сферы (исключить перекрытие потенциалов, отвечающих разным атомам), а в пространстве между сферами положить потенциал равным некоторой константе. Коэффициенты перед этими функциями подбираются таким образом, чтобы волновая функция на поверхности сферы, переходя в плоскую волну, не испытывала скачка. Присоединённые плоские волны, как и ортогонализованные, обеспечивают достаточно быструю сходимость ряда и образуют полную систему функций. 5.1.2. Орбитали Слэтера Базисы гауссовых и слэтеровских орбиталей применяется, в основном, для расчета конечных тел: кластеров и отдельных молекул. Слэтеровские орбитали являются точным решением для водородоподобного атома и наиболее точны, если расчёт начинается с изолированных атомов. Однако некоторые интегралы кулоновского типа не могут быть решены аналитически при использовании этого набора. Наиболее распространены в квантово-химических расчётах простейшие орбитали слэтеровского типа — атомные орбитали Слэтера—Зенера: s 2ζ (2ζ)n rn−1 e−ζr Ylm (θ, φ), χSZ nml (ζ, r) = (2n)! где  cos(mφ), если m > 0, 2l + l(l − |m|!) |m| P (cos θ) Ylm (θ, φ) = sin(mφ), если m < 0, 2πδm (l − |m|!) l — нормированные вещественные сферические функции, |m| Pl (cos θ) |m| 1 dl+|m| 2 2 = l l − (cos θ) (cos2 θ − 1) l+|m| 2 l! (d cos θ) 45

— присоединенные полиномы Лежандра,  2, если m = 0, δm = 1, иначе. (5.2) Параметрами орбиталей являются квантовые числа n, l и m, показатели экспоненты ζ и координаты точек центрирования орбиталей. Обычно для базисных функций используются квантовые числа, соответствующие занятым и ближайшим возбужденным состояниям электрона в атоме. В целях сокращения количества молекулярных интегралов можно включать в базисный набор только орбитали с главными квантовыми числами, не превосходящими главного квантового числа n для заполненных оболочек в атоме (минимальный базис). В случае использования полуэмпирических методов в базис включаются только валентные орбитали (валентное приближение) [33]. Для определения ζ Слэтер в 1930 году предложил набор эмпирических правил, основанных на стремлении к наилучшему воспроизведению первых потенциалов ионизации атомов. Согласно этим правилам Z −S , ζ= n∗ где Z — заряд ядра, n∗ — эффективное (возможно, дробное) главное квантовое число, S — константа экранирования. Константа экранирования S равна сумме приписываемых каждой группе электронов в атоме (1s, 2s, 2p,. . . ) констант sm : S= X sm , m sm   0        0,35     0,3 =   0,85       1     1 для электронов внешних групп, для каждого электрона в группе, кроме рассматриваемого, для электрона в 1s группе, кроме рассматриваемого, для каждого из s- и p-электронов с n = n0 − 1, для каждого электрона следующих внутренних оболочек, для электронов внутренних d- и f - групп (здесь n — главное квантовое число рассматриваемого электрона). 0 46

Радиальные функции слэтеровского типа недостаточно точно описывают поведение атомных орбиталей Хартри—Фока на небольших расстояниях от ядра. Этот недостаток устраняется аппроксимацией каждой из этих атомных орбиталей по крайней мере двумя слэтеровскими функциями с разными орбитальными экспонентами — дубль-зета-функциями: n∗−1 −ζ1 r n∗−1 −ζ2 r χSZ = N r e + r e Ylm (θ, φ). nlm Квантово-химические расчеты с использованием слэтеровских функций в качестве базисных орбиталей являются достаточно сложными и трудоемкими. Гауссовы орбитали частично решают эту проблему, но иногда при этом разложении также необходимо достаточно большое число членов в ряду. 5.1.3. Орбитали Гаусса В 1950 году С. Ф. Бойс предложил использовать базисные функции гауссова типа, имеющие в декартовых координатах вид 2 lx ly lz −αr χG , nlm (r, a, l) = x y z e где lx , ly , lz — целые положительные числа, l = lx + ly + lz . В сферических координатах они имеют вид r χG nlm = n+1/2 2 (4α) π (2n − 1)!! !1/2 2 rn−1 e−αr Ylm (θ, φ), 2 где Ylm — сферические гармоники, rn−1 e−αr — радиальные функции типа Гаусса. Их основное отличие от слэтеровских функций заключается в квадратичной зависимости от аргумента r экспоненты. Это позволяет представлять произведения гауссовых функций (центрированных в разных точках), встречающиеся в многоцентровых интегралах, в виде линейной комбинации гауссовых функций, центрированных в общем случае в другой точке пространства. Именно для этого класса функций оказывается возможным 47

аналитическое вычисление многоцентровых интегралов, поскольку произведение двух любых функций типа Гаусса с различными точками отсчёта также является гауссовой функцией: 2 2 2 e−am |r−rm | e−at |r−rt | = smt e−amt |r−rmt | , amt = am +at , rmt = (am rm + at rt ) . amt Множитель smt быстро убывает с увеличением расстояния rm − rt , благодаря чему значительное число произведений базисных функций вносит пренебрежимо малый вклад и может не вычисляться. Это упрощает вычисления, поэтому одну слэтеровскую атомную орбиталь, как правило, аппроксимируют несколькими гауссовыми функциями. Наиболее простым типом базисных наборов является ОСТ-nГФ (ST O-nG), где каждая атомная орбиталь состоит из суммы n (обычно от двух до шести) функций типа Гаусса, причём коэффициенты подобраны с целью приближения поведения линейных комбинаций к орбиталям типа Слэтера. Проведение тестовых расчетов с использованием этих базисных наборов показало, что при n ≥ 3 результаты расчётов очень схожи. Поэтому наиболее широкое распространение получил минимальный базисный набор ОСТ-3ГФ (ST O-3G). Минимальные базисные наборы включают только атомные орбитали, которые необходимы для размещения электронов нейтрального атома. Сферическая симметрия атомов и пространственная инвариантность молекул требуют включение всех трёх типов np-орбиталей при появлении хотя бы одного p-электрона. Малый размер и простота базиса ОСТ-3ГФ — его основное преимущество. Однако необходимо отметить недостатки: неудовлетворительные результаты расчётов соединений электроположительных элементов третьего периода, переоценка стабильности малых циклов и π-акцепторной способности электроположительных элементов второго периода. Эти проблемы устранимы использованием более широких базисов: валентно-расщепленных и биэкспоненциальных. В них атомные орбитали представлены двумя частями: внутренней и внешней (при построении молекулярных орбиталей коэффициенты Cij можно варьировать независимо). 48

В валентно-расщепленных базисных наборах на внутреннею и внешнюю составляющие разделены только валентные орбитали. Поэтому схема записывается как m-npГФ (m-npG), где m — число гауссовых функций, заменяющих каждую внутреннюю атомную орбиталь, n и p — число гауссовых функций с разными значениями экспонент, аппроксимирующих каждую валентную атомную орбиталь. Наиболее часто применяемым является базис 6-31ГФ (6-31G), который фигурирует и в данной работе. В биэкспоненциальных базисах расщеплены и валентные, и внутренние орбитали. Для улучшения описания молекулярных орбиталей на больших расстояниях от ядра часто используют базисные наборы с включением поляризационных функций. Поляризационные функции — дополнительные волновые функции с l + 1 орбитальным квантовым числом последних заполненных атомных орбиталей. Одними из наиболее используемых поляризационных базисных наборов являются 6-31ГФ∗ и 6-31ГФ∗∗ , где звёздочка обозначает добавление поляризационных d-функций к p-элементам (или f -функций к d-элементам). Вторая звёздочка обозначает добавление поляризационных p-функций к 1s-орбиталям атомов водорода. Для расчётов молекул, требующих более точного описания несвязывающих электронных пар, в базисные наборы вводятся специальные диффузные s- и p-функции со значениями экспонент от 0,1 до 0,01. Их включение обозначается знаком «+», например, 3-21+ГФ. Выбор базиса, подходящего для решения задачи, определяется стремлением получить более точное решение, с одной стороны, и ограничениями, связанными с ресурсами ЭВМ, с другой. Поэтому на практике обычно полная оптимизация выполняется с использованием небольших базисов, после чего в более широких базисах проводятся расчёты для одной геометрической конфигурации, и устанавливаются поправки, связанные с учётом электронной корреляции. 49

5.2. Численное интегрирование выражений для вычисления обменно-корреляционной энергии Так как обменно-корреляционный функционал сложным образом зависит от переменной интегрирования (через электронную плотность), найти этот интеграл можно только численно. Успешное применение метода функционала электронной плотности стало возможным благодаря изобретению подходящих схем трёхмерного численного интегрирования. Как правило, пространственная сетка строится с количеством узлов от 1000 до 30000 на каждый атом [34]. В программе Gaussian сетка интегрирования по умолчанию состоит из 75 радиальных оболочек с 302 угловыми точками на каждую оболочку. Чем больше сетка, тем меньше шаг интегрирования и точнее результат численного интегрирования, поэтому значение энергии системы зависит и от выбора сетки. Рассмотрим некоторый интеграл Z I = f (r)dr, (5.3) где, например, f (r) = εhyb xc . Его вычисление затруднено быстрым изменением f (r) в окрестности каждого ядра или отсутствием непрерывности в точке его расположения. Одним из способов решения этой задачи является разбиение пространства молекулы на неперекрывающиеся области, в каждой из которых находится не более одного ядра. Более простой подход состоит в использовании атомных весовых функций ωk (r), близких к единице в окрестности нуля, стремящихся к нулю по мере приближения к другому ядру и удовлетворяющих условию X ωK (r) = 1. K Такое представление позволяет разбить «многоцентровый» интеграл (5.3) на «одноцентровые»: I= X K 50 IK ,

Z IK = Z ωK (r)f (r)dr = fK (r)dr, где fK (r) локализована около соответствующего K-ого ядра. Каждый одноцентровый интеграл с центром в точке расположения соответствующего ядра удобно берётся в сферических координатах путём раздельного численного интегрирования по радиальной r и угловым s переменным: Z ∞ r2 IK = Z fK (RK + rs)dsdr. |s|=1 0 Атомные весовые функции ищут в виде PK (r) , ωK (r) = P P (r) J J где в качестве PK (r) выступают различные локализованные функции, например, так называемые «ячеечные» |r − RK | − |r − RJ | PK (r) = ΠJ6=K s3 , |RK − RJ | 1 1 s3 (t) = (1 − p(p(p(t)))), p(t) = (3t − t3 ), 2 2 многочлен s3 — «сглаженный» аналог ступенчатой функции  1, −1 ≤ t ≤ 0, s(t) = 0, 0 < t ≤ 1. Вычислительные затраты на построение обменно-корреляционного функционала в базисе функций типа Гаусса с помощью численного интегрирования за счёт вычисления в каждой точке сетки лишь значений базисных функций ближайшего окружения приводят к линейной зависимости от размера задачи [34]. Перечислим также основные источники погрешностей квантовохимических расчётов. Среди них необходимо назвать: • предположение об отсутствии релятивистского эффекта — малости скорости движения электронов по сравнению со скоростью света 51

(справедливо с хорошей степенью точности для атома водорода и элементов II периода, но существенно для внутренних электронов элементов IV–VI периодов), • описание движения электронов в поле покоящихся ядер (приближение Борна—Оппенгеймера), • несовершенство базисного набора, • приближённое описание коррелированного движения электронов, • приближённое вычисление интегралов и ошибки округления. 5.3. Запуск задач в Gaussian и распределённые вычисления После выбора начального приближения для волновой функции (в базисе 6-31g) и конкретной реализации метода функционала электронной плотности (с приближением для обменно-корреляционного функционала B3LYP) можно приступать к расчётам. Вычисления были выполнены на кластере высокопроизводительных вычислений факультета прикладной математики — процессов управления и на виртуальной машине ресурсного центра Научного парка СПбГУ «Вычислительный центр». Кластер является разновидностью распределённой (параллельной) системы, состоящей из нескольких соединённых между собой компьютеров, использующихся как единый ресурс. Удалённое управление было осуществлено посредством сетевого протокола SSH (англ. secure shell — безопасная оболочка). Он позволяет безопасно передавать другие протоколы в незащищённой среде. В отличие от других протоколов (telnet, rlogin и др.) шифрует весь трафик, включая пароли, а также позволяет предварительно сжимать данные (по запросу). Кроме работы в командной строке возможен удалённый запуск программ, в том числе графических. Для работы по SSH нужен SSH-клиент и SSH-сервер. Клиент используется для входа на удалённый компьютер и выполнения команд. Сервер после связи с ним проводит аутентификацию и начинает обслуживание. 52

Для соединения клиента с сервером оба должны создать пары открытых и закрытых ключей и обменяться открытыми. Как правило, также используется пароль. Подключение осуществлялось посредством клиента PuTTY. Для работы с графическими приложениями удалённого компьютера использовался сервер Xming. Моделирование объектов исследования проводилось с помощью программы Chemcraft, расчёты были выполнены в квантово-химическом пакете Gaussian. Последний предлагает реализации большого количества методов квантовой и вычислительной химии. Программа Gaussian была выпущена в 1970 г. группой под руководством нобелевского лауреата по химии (1998 г.) Дж. Поплом; на сегодняшний день считается наиболее распространённым и мощным инструментом моделирования. Использовалась реализация Gaussian 09 и графическое приложение GaussView 5. Опишем схему взаимодействия с программой Gaussian. На вход подаётся файл с расширением .com, в котором обязательно указывается расчётный метод, выбранный базис функций (можно задавать самостоятельно или пользоваться готовыми наборами) и начальные координаты системы с указанием типа атома (или атомным номером), зарядом и мультиплетностью. Затем с помощью ключевых слов указывается тип задачи, например, нахождение оптимальной геометрии (ключевое слово opt), расчёт частот колебаний (freq) и др. Можно также варьировать количество выделяемой под задачу памяти. Пример входного файла с декартовыми координатами: %mem=10GB // выделеление памяти %chk=c30.chk // создание chk-файла #p rb3lyp/6-31g opt // метод, базис, тип задачи C(3,0) // заголовок файла 0 1 // заряд и мультиплетность 53

C 1.216023 0.000000 -11.516243 C 1.216023 0.000000 -12.937243 C 0.608011 1.053107 -9.360036 С ... //координаты атомов x,y,z Результатом работы программы является log-файл. Его содержимое отличается в зависимости от задачи, но содержит информацию о всех операциях. Если это было указано в исходном .com файле, то будет создан и chk-файл — файл с сохранением контрольной точки. Приведём в пример некоторые типы задач и соответствующие им ключевые слова. Opt Ключевое слово, запускающее процесс оптимизации геометрии, который циклически продолжается, пока не будет найдена стационарная точка поверхности потенциальной энергии (локальный минимум энергии). По умолчанию оптимизация проводится во внутренних координатах. Математически это осуществляется следующим образом. В методе функционала электронной плотности после перехода от стационарного уравнения Шрёдингера к задаче на экстремум функционала, имеющего смысл оператора энергии, получается система алгебраических уравнений, именуемых уравнениями Кона—Шэма. Пробная волновая функция, формируемая с помощью заданного пользователем базиса, подставляется в эту систему уравнений и решается заново, пока не будут достигнуты критерии сходимости. Критерием сходимости является обращение в нуль вектора ∇E: ∂E ∂E ∂E ∂E ∂E ∂E ∇E = , , , . . . , 1 , 2 , 3 = 0, ∂x11 ∂x21 ∂x31 ∂xK ∂xK ∂xK где K — число атомов в структуре, x1 , x2 , x3 — декартовы координаты. С точки зрения приближения Борна—Оппенгеймера E(R) — потенциальная энергия ядер, а -∇E представляет обобщённую силу, поэтому нулевыми должны быть и силы, действующие на каждый атом. Gaussian реализует следующие критерии следующие: • ограничение на максимальное значение компоненты силы 54

− ∂E ! < 0,00045 a.e., ∂xji • ограничение на максимальное смещение координаты по сравнению со max i=1,...,K; j=1,2,3 значением на предыдущем шаге ∆xj ≤ 0,0018 Å, • среднеквадратичное отклонение всех сил от среднего значения < 0,0003 а. е., • среднеквадратичное отклонение всех смещений от среднего смещения ≤ 0,0012 Å. Иногда достичь сходимости не удаётся: значения энергии колеблются около предположительно стационарной точки. В этом случае можно попытаться воспользоваться улучшенными критериями сходимости, например, квадратичным. Управление процедурой сходимости осуществляется с помощью ключевого слова SCF. Другой способ — принудительное увеличение количества итераций или ухудшение критерия сходимости. Также в некоторых случаях удобно проводить расчёты на неоднородной сетке, когда известно, что на одном участке системы требуется высокая точность. Это осуществляется с помощью технологии ONIOM расчётов, предложенной японским профессором К. Морокумой. Однако при нарушении сходимости необходимо помнить, что эта процедура позволяет найти только локальный минимум энергии. Cubegen Генерирует cube-файлы, которые позволяют визуализировать данные (электронную плотность, электростатический потенциал и др.). Сube-файл представляет собой файл с трёхмерной равномерной сеткой. Formchk Утилита, с помощью которой осуществляется преобразование сформированного программой chk-файла в текстовый fchk-файл, имеющий те же данные в удобном представлении. 55

6. Результаты численного эксперимента 6.1. Результаты моделирования структур Первая часть моделирования состояла в нахождении координатных представлений структур и получении оптимальных в смысле полной энергии геометрий. В результате оптимизации геометрии три модели с меньшим количеством атомов C72 (3, 0) (рис. 6.1), C72 (3, 3) (рис. 6.2), Si36 C36 (3, 3) (рис. 6.3) образовали замкнутые структуры. Рис. 6.1: Оптимизированная структура C72 (3, 0) Рис. 6.2: Оптимизированная структура C72 (3, 3) Рис. 6.3: Оптимизированная структура Si36 C36 (3, 3) Наконечник нанотрубок C120 (5, 5) и C240 (5, 5) (рис. 6.4) искривился, но не замкнулся. Этот результат косвенно подтверждает тот факт, что 56

нанотрубки, диаметр которых больше или равен характерному для (5, 5) диаметру, могут быть «закрыты» частями фуллерена [35]. Рис. 6.4: Оптимизированная структура C240 (5, 5) Характерные геометрические параметры структур приведены в таб. 6.1. Таблица 6.1: Длина и диаметр объектов изучения Модель C72 (3, 0) C72 (3, 3) Si36 C36 (3, 3) C120 (5, 5) C240 (5, 5) Диаметр, Å 2,35 4,07 4,07 6,78 6,78 Длина, Å 23,28 14,76 14,76 14,76 29,51 Процедура самосогласования в среднем занимала 10 циклов. Другими словами, количество итераций в методе функционала электронной плотности, и, соответственно, количество улучшений пробной волновой функции, было в среднем равно 10. Частично этого сравнительно небольшого значения удалось добиться благодаря предварительной оптимизации структур. Но даже в этом случае время, затраченное на решение уравнений, доходило до нескольких суток. На рис. 6.5 приведён пример графика сходимости энергии от числа итераций, вид которого не свидетельствует против корректности решения. 57

Рис. 6.5: График сходимости энергии для C72 (3, 3) Вторым этапом работы было моделирование выбранных структур во внешнем электрическом поле. Влияние этого поля на объекты изучения оценивалось по изменениям распределения электронной плотности в структурах, плотности тока эмиссии, полных электронных энергий, дипольных моментов и ширины запрещённых зон. 6.2. Распределение электронной плотности Следствием воздействия внешнего электрического поля на структуры является, в первую очередь, перераспределение электронной плотности. Ожидаемым результатом является симметричное распределение в отсутствии внешнего поля и асимметричное — в противном случае. 58

Наглядно продемонстрировать этот факт позволяют заряды Малликена. Они отражают картину распределения электронной плотности путём равного разделения последней между двумя атомными орбиталями химически связанных атомов. Эти заряды являются условными и определяются как разность между зарядом ядра и полной атомной электронной заселенностью атома, то есть в виде qA = zA − X Pij Sij , i,j∈A где zA — заряд ядра, Pij = P k cki ckj — матричные элементы зарядов и порядков связей, Sij — матричные элементы интегралов перекрывания базисных функций. Из этого определения следует, что значения зарядов Малликена зависят от выбора базисного набора. При использовании больших наборов базисных функций или добавлении диффузных функций они могут давать физически необоснованные значения. Базис 6-31g, использованный в работе, является минимальным и не включает диффузные функции. Схема отнесения электронной плотности к атомам, предложенная американским химиком Р. С. Малликеном, при некоторых недостатках позволяет судить о распределении электронной плотности и влиянии внешнего электрического поля на структуры. На рис. 6.6–6.8 приведены картины распределения зарядов Малликена в отсутствии внешнего электрического поля (верхняя часть каждого изображения с симметричным относительно центра системы распределением зарядов) и во внешнем электрическом поле напряжённости 0,25 B/Å (нижняя часть с асимметричным распределением). Цвета и приведённые значения отвечают минимальному (синий) и максимальному (красный) значению заряда в системе. Во всех случаях наблюдается смещение электронной плотности в направлении вектора напряжённости поля. 59

Рис. 6.6: Распределение зарядов Малликена для C72 (3, 0) вне и во внешнем электрическом поле напряжённости 0,25 B/Å Рис. 6.7: Распределение зарядов Малликена для Si36 C36 (3, 3) вне и во внешнем электрическом поле напряжённости 0,25 B/Å 60

Рис. 6.8: Распределение зарядов Малликена для C240 (5, 5) вне и во внешнем электрическом поле напряжённости 0,25 B/Å 6.3. Плотность эмиссионного тока Решение стационарного уравнения Шрёдингера подразумевает нахождение полной энергии основного состояния и соответствующего значения волновой функции. Исходя из предположения о способе формирования пробной волновой функции её значение можно найти, зная коэффициенты в разложении (5.1). Отсюда на основании формулы (4.3) можно оценить плотность эмиссионного тока. Эти оценки приведены в таб. 6.2. Таблица 6.2: Значения плотности эмиссионного тока нанотрубок углерода и карбида кремния вне поля (j) и в поле напряжённости 0,25 B/Å (jext ), полученные из численного эксперимента Модель C72 (3, 0) C72 (3, 3) Si36 C36 (3, 3) C120 (5, 5) C240 (5, 5) j, пА/Å 2,61 0,40 20,55 0 4,41 jext , пА/Å 29,79 53,84 211,78 87,82 213,84 61

Очевидно, плотность тока эмиссии в отсутствии внешнего электрического поля должна быть близка к нулевой. Наибольшее отклонение от нуля в этом случае имеет нанотрубка карбида кремния Si36 C36 (3, 3), что может быть связано с неточностью описания соединения атомов двух разных типов — кремния и углерода. Кроме того, базис 6-31g, выбранный в качестве формирующего пробную волновую функцию, наиболее точно позволяет описывать поведение атомов II периода, к которым относится углерод. Кремний же является элементом III периода. Однако все расчёты были проведены с использованием базиса 6-31g, так как сопоставление результатов возможно только в рамках одного метода и его реализации. Характер изменения величин в целом оправдывает ожидаемый: в отсутствии внешнего поля это величины, близкие к нулю, которые увеличиваются с возрастанием напряжённости приложенного электрического поля (рис. 6.9). 250 C72(3,0) C72(3,0) Si36C36(3,3) 200 C120(5,5) C240(5,5) j, pA/A' 150 100 50 0 0,00 0,05 0,10 0,15 0,20 0,25 F, B/A' Рис. 6.9: График зависимости плотности эмиссионного тока j от напряжённости электрического поля F 62

Для всех углеродных структур во внешнем электрическом поле наблюдается некоторая зависимость плотности тока эмиссии от количества атомов. Однако эта зависимость не линейна: две структуры одной конфигурации (5, 5), длины которых отличаются в два раза (C120 и C240 ), испускают не в отношении 1:2. Тип симметрии структур, вероятно, также может влиять на эмиссионные характеристики. Углеродные структуры различных типов симметрий — зигзагообразной (3, 0) и кресельной (3, 3) — излучают не одинаково, хотя содержат равное количество атомов. Рассмотренная нанотрубка карбида кремния короче и у́же структуры C240 (5, 5), однако по результатам численного эксперимента имеет близкую эмиссионную способность. 6.4. Полные энергии и дипольные моменты систем О влиянии электрического поля на структуры можно судить по изменению полной электронной энергии и дипольного момента систем. Напряжённость поля в 0,25 B/Å в атомных единицах соответствует 0,005 а.е.н. Как видно из таб. 6.3, наложение поля предсказуемо способствует уменьшению полных энергий систем. Для углеродных трубок можно проследить линейную зависимость полной энергии от количества атомов: полные энергии убывают с возрастанием количества атомов. Этот результат также свидетельствует о корректности моделирования. Таблица 6.3: Значения полных энергий систем в отсутствии поля (E) и во внешнем электрическом поле напряжённости 0,005 а.е.н. (Eext ), полученные из численного эксперимента Модель C72 (3, 0) C72 (3, 3) Si36 C36 (3, 3) C120 (5, 5) C240 (5, 5) −2739,81 −2741,57 −11791,40 −4570,67 −9143,18 Eext , а.е. −2739,73 −2741,59 −11791,42 −4570,69 −9143,29 E, а.е. Результирующие значения дипольных моментов для каждого случая приведены в таб. 6.4. 63

Таблица 6.4: Значения дипольных моментов систем в отсутствии поля (p) и во внешнем электрическом поле напряжённости 0,005 а.е.н. (pext ), полученные из численного эксперимента Модель C72 (3, 0) C72 (3, 3) Si36 C36 (3, 3) C120 (5, 5) C240 (5, 5) p, Д 0 0 0,008 0 0,0002 pext , Д 65,49 19,69 31,82 27, 83 112 Известно, что дипольный момент неполярных систем в отсутствии внешнего электрического поля должен быть равен нулю. Это с высокой точностью выполняется для всех рассматриваемых углеродных нанотрубок. Для карбидокремниевой нанотрубки эта величина незначительно превышает погрешность и может быть связана с присутствием атомов кремния. В поле напряжённости 0,005 а.е.н. наибольший дипольный момент наблюдается у трубки с наибольшим количеством атомов, однако среди трёх моделей C72 (3, 0), C72 (3, 3), Si36 C36 (3, 3) с одинаковым количеством атомов наиболее подвержена влиянию поля трубка наименьшего диаметра C72 (3, 0). При одинаковых конфигурациях структур Si36 C36 (3, 3) и C72 (3, 3) дипольный момент выше у первой. 6.5. Анализ зонной структуры нанотрубок В рамках метода молекулярных орбиталей был проведён анализ зонной структуры моделей. Значения ширины запрещённой зоны в электронвольтах указаны в таб. 6.5. Полученные для Si36 C36 (3,3) значения ширины запрещённой зоны очень близки к экспериментальным [36] и согласуются с фактом, что карбид кремния — широкозонный полупроводник. Результирующие значения ширины запрещённой зоны углеродных нанотрубок также принадлежат интервалу, характерному для этого материала, однако нуждаются в комментариях. В соответствии с приведённым в п. 3.2.1 тезисом о связи зонной 64

Таблица 6.5: Значения ширины запрещённой зоны нанотрубок углерода и карбида кремния вне поля (∆E) и в поле напряжённости 0,25 B/Å (∆Eext ), полученные из численного эксперимента Модель C72 (3, 0) C72 (3, 3) Si36 C36 (3, 3) C120 (5, 5) C240 (5, 5) ∆E, эВ 0,95 1,37 2,37 0,86 1,01 ∆Eext , эВ 0,90 1,26 2,34 0,80 0,90 структуры с конфигурацией нанотрубок, рассмотренные модели должны обладать металлическим типом проводимости. Результаты проведённого численного эксперимента свидетельствуют, что данные модели прямых углеродных трубок являются полупроводниками. Ниже приведены диаграммы орбитальных энергий для карбидокремниевой (рис. 6.10) и углеродных (рис. 6.11) нанотрубок, отражающие их зонную структуру. Рис. 6.10: Зонная структура Si36 C36 (3, 3) в основном состоянии 65

C72 (3, 0) C72 (3, 3) C120 (5, 5) C240 (5, 5) Рис. 6.11: Зонная структура моделей углеродных нанотрубок в основном состоянии 66

7. Заключение Целью работы было теоретическое исследование влияния внешнего электрического поля на различные структуры открытых одностенных нанотрубок углерода и карбида кремния. Полученные результаты необходимы для создания математической и компьютерной модели полевой электронной эмиссии на нано-уровне. Обзор литературы позволил убедиться в важности изучения нанотрубок как возможного будущего электроники. Результаты подтвердили факт, что структуры из карбида кремния аналогичных с углеродными конфигураций более подвержены влиянию сильного электрического поля. При условии бо́льшей чем у углеродных структур стабильности и долговечности нанотрубки карбида кремния можно считать не менее перспективным материалом. На эмиссионные характеристики влияет диаметр нанотрубок, однако можно предположить, что бо́льшее значение с этой точки зрения имеет общее количество атомов в структуре. Мы подробно рассмотрели и обосновали применение метода функционала электронной плотности для моделирования поведения нанотрубок. Были освещены альтернативные формы моделирования и отмечены преимущества выбранного метода как обеспечивающего наиболее корректное описание реальных физических явлений на квантовом уровне. Для практической реализации задачи были изучены возможности квантово-химического пакета Gaussian 09 и основы теории распределённых вычислительных систем, освоена технология удалённого подключения к кластеру высокопроизводительных вычислений. Были получены координатные представления систем, оптимизирована по энергии геометрия объектов, выбраны начальное приближение для волновой функции (в базисе 631g) и конкретная реализация метода функционала электронной плотности (B3LYP). В результате были получены и оценены эмиссионные характеристики нанотрубок углерода и карбида кремния различных конфигураций, сделан вывод об адекватности моделирования. В заключение необходимо сделать ряд важных уточнений. Нами бы67

ли рассмотрены идеальные структуры нанотрубок. В реальности, кроме проблемы получения однородного массива нанотрубок в достаточном количестве, имеют место множественные дефекты структур, такие как вакансии или локтевые соединения; трубки различных конфигураций могут быть соединены в одну, в многослойных нанотрубках могут отсутствовать части слоёв. Всё это осложняет создание общей модели поведения нанотрубок. Кроме того, в силу своего масштаба они не могут существовать без подложки, граница раздела с которой также может влиять на процесс эмиссии. Исследования были проведены с использованием оборудования ресурсного центра Научного парка СПбГУ «Вычислительный центр» и кластера высокопроизводительных вычислений факультета ПМ—ПУ. 68

Список литературы 1. Егоров Н. В., Шешин Е. П. Автоэлектронная эмиссия. Принципы и приборы. М.: Интеллект, 2011. 704 с. 2. Fan J., Chu P. K. Silicon Carbide Nanostructures: Fabrication, Structure, and Properties. Berlin: Springer, 2014. 330 p. 3. Ilyin V. A., Luchinin V. V. et al. Superfast drift step recovery diodes (DSRDs) and vacuum field emission diodes based on 4H-SiC // Materials Science Forum. 2013. Vol. 740–742. P. 1010–1013. 4. Gogotsi Yu., Presser V. Carbon Nanomaterials. Second Edition. Florida: CRC Press, 2014. 512 p. 5. Дьячков П. Н. Углеродные нанотрубки: строение, свойства, применения. М.: БИНОМ. Лаборатория знаний, 2006. 293 с. 6. Гусев А. И. Наноматериалы, наноструктуры, нанотехнологии. М.: Физматлит, 2007. 416 с. 7. Елинсон М. И., Васильев Г. Ф. Автоэлектронная эмиссия. М.: Физматлит, 1958. 272 с. 8. Li Z. Density functional theory for field emission from carbon nanostructures // Ultramicroscopy. 2015. Vol. 159. No 2. P. 162–172. 9. Попов А. M. Вычислительные нанотехнологии. М.: Издательский отдел факультета ВМК МГУ им. М. В. Ломоносова, 2009. 280 c. 10. Ding F., Jiao K., Wu M., Yakobson B. I. Pseudoclimb and dislocation dynamics in superplastic nanotubes // Phys. Rev. 2007. Vol. 98. No 7. P. 1435–1438c. 11. Фёдоров А. С., Сорокин П. Б., Аврамов П. В., Овчинников С. Г. Моделирование свойств, электронной структуры ряда углеродных и неуглеродных нанокластеров и их взаимодействия с лёгкими элементами. Новосибирск: СО РАН, 2006. 435 с. 12. Novoselov K. S. Electric field effect in atomically thin carbon films // Science. 2004. Vol. 306. P. 34–39. 13. Dresselhaus M. S., Dresselhaus G., Saito R. Physics of carbon nanotubes // Carbon. 1995. Vol. 33. No 7. P. 883–891. 69

14. Мурзашев А. И., Шадрин Е. О. Энергетический спектр и оптические свойства бесконечных углеродных нанотрубок в модели Хаббарда // Физика твердого тела. 2012. Т. 54. № 12. С. 2359–2365. 15. Peng J., Li Z., He C., Chen G., Wang W., Deng Sh., Xu N., Zheng X., Chen G., Edgcombe C. J., Forbes R. G. The roles of apex dipoles and field penetration in the physics of charged, field emitting, single-walled carbon nanotubes // J. Appl. Phys. 2008. Vol. 104. P. 289–302. 16. Zheng X., Chen G., Li Z., Deng Sh., Xu N. Quantum-mechanical investigation of field-emission mechanism of a micrometer-long singleWalled carbon nanotube // Phys. Rev. 2004. Vol. 92. P. 125– 129. 17. Chen G., Wang W., Peng J., He C., Deng Sh., Xu N., Li Z. Atomic decoration for improving the efficiency of field electron emission of carbon nanotubes // Phys. Rev. 2007. Vol. 76. P. 173–179. 18. Choyke W. J., Pensl G. Physical properties of silicon carbide // MRSBulletin. 1997. P. 25–29. 19. Svetlichnyi A. M., Spiridonov O. B., Volkov E. Y., Linets L. G., Grigoriev M. N. The evaluation of the characteristics of field emission nanostructures based on Si and SiC // Izvestiya SFedU. Engineering sciences. 2011. P. 27–34. 20. Охрименко О. Б., Конакова Р. В., Светличный А. М., Спиридонов О. Б., Волков Е. Ю. Оценка автоэмиссионных свойств наноструктур на основе карбида кремния и графена // Наносистемы, наноматериалы, нанотехнологии. 2012. Т. 10. № 2. С. 335–342. 21. Овсянникова Л. И., Покропивный В. В., Бекенев В. Л. Электронная структура кристаллообразующих фуллеренов, фулсиценов и кристаллов из них — фулсиценитов // Физика твёрдого тела. 2009. Т. 51. № 10. C. 2070–2077. 22. Левич В. Г. Курс теоретической физики. Том II. М.: Наука, 1971. 912 с. 23. Глинка Н. Л. Общая химия. М.: Интеграл-пресс, 2005. 728 с. 24. Фок В. А. Начала квантовой механики. М.: Наука, 1976. 376 с. 25. Ландау Л. Д., Лифшиц Е. М. Квантовая механика. М.: Наука, 1963. 70

767 с. 26. Hohenberg P., Kohn W. Inhomogeneous electron gas // Phys. Rev. B. 1964. Vol. 136. P. 1182–1196. 27. Хартри Д. Расчёты атомных структур. М.: ИИЛ, 1960. 256 с. 28. Зеленцов С. В. Экспериментальные и теоретические методы изучения возбужденных состояний. Н. Новгород: Изд-во Нижегородского гос. ун-та им. Н. И. Лобачевского, 2007. 78 с. 29. Kohn W., Sham L. J. Self-consistent field equations including excange and correlation effects // Phys. Rev. A. 1965. Vol. 140. P. 1133–1138. 30. Lee C., Yang W., Parr R. G. Development of the Colle-Salvetti correlation energy formula into a functional of the electron density // Phys. Rev. B. 1988. Vol. 37. P. 785–789. 31. Becke A. D. Density-functional exchange-energy approximation with correct asymptotic behavior // Phys. Rev. A. 1988. Vol. 38. P. 3098–3100. 32. Мартинсон Л. К., Смирнов Е. В. Квантовая физика. М.: Изд-во МГТУ им. Н. Э. Баумана, 2006. 528 с. 33. Крайнов В. П. Лекции по микроскопической теории атомного ядра. М.: Атомиздат, 1973. 224 с. 34. Лайков Д. Н. Развитие экономного подхода к расчёту молекул методом функционала электронной плотности и его применение к решению сложных химических задач. М.: Изд-во МГУ им. М. В. Ломоносова, 2000. 103 с. 35. Saito R., Fujita M., Dresselhaus M. S., Dresselhaus G. Electronic structure of graphene tubules based on C60 // Phys. Rev. B. 1992. Vol. 46. No 3. P. 1804. 36. Елецкий А. В. Транспортные свойства углеродных нанотрубок // Успехи физических наук. 2009. Т. 179. № 3. С. 209–224. 71

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв