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

Чепелев И.Г.

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

22.04.01 Материаловедение и технологии материалов

Общие и комплексные проблемы технических и прикладных наук и отраслей народного хозяйства

Дипломы

Вуз: Белгородский государственный университет - национальный исследовательский университет (НИУ «БелГУ»)

ID: 5c1a53417966e104f6f852a8

UUID: e8e8f4d0-e5c6-0136-fcf0-525400005860

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

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

Просмотры: 9

Чепелев И.Г.

Источник: Федеральное государственное автономное образовательное учреждение высшего образования «Белгородский государственный национальный исследовательский университет»

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

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

Скачать - 1,5 МБ

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

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ «БЕЛГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ» (НИУ «БелГУ») ИНСТИТУТ ИНЖЕНЕРНЫХ ТЕХНОЛОГИЙ И ЕСТЕСТВЕННЫХ НАУК КАФЕДРА МАТЕРИАЛОВЕДНИЯ И НАНОТЕХНОЛОГИЙ ВЛИЯНИЕ ТОЧНОСТИ ЗАДАНИЯ ТРЕХЧАСТИЧНЫХ ВЗАИМОДЕЙСТВИЙ НА ПРОГНОЗ ХАРАКТЕРИСТИК СТРУКТУР УГЛЕРОДА Выпускная квалификационная работа обучающегося по направлению подготовки Материаловедение и технологии материалов 22.04.01 очной формы обучения, группы 07001639 Чепелева Игоря Геннадьевича Научный руководитель доктор физико-математических наук Липницкий А.Г. БЕЛГОРОД 2018

Содержание Введение 4 1 Обзор литературы 7 1.1 Теория функционала электронной плотности . . . . . . . . . 7 1.1.1 Псевдопотенциалы . . . . . . . . . . . . . . . . . . . . 15 1.1.2 Расчеты электронной структуры . . . . . . . . . . . . 17 1.2 Теория динамики решетки . . . . . . . . . . . . . . . . . . . 21 1.2.1 Модель Борна-вон Кармана . . . . . . . . . . . . . . . 21 1.2.2 Метод сверхячеек . . . . . . . . . . . . . . . . . . . . . 24 1.2.3 Вклад тепловых колебаний в термодинамические свойства твердых тел . . . . . . . . . . . . . . . . . . . . . 26 1.3 Молекулярная динамика . . . . . . . . . . . . . . . . . . . . . 27 1.3.1 Термодинамические ансамбли . . . . . . . . . . . . . . 32 1.3.2 Молекулярная статика . . . . . . . . . . . . . . . . . . 33 1.4 История развитие межатомных потенциалов углеродных систем . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.5 Межатомные потенциалы в рамках нового подхода . . . . . 41 1.6 Процедура оптимизации параметров потенциалов межатомных взаимодействий . . . . . . . . . . . . . . . . . . . . . . . 2 Материал и методика исследования 2.1 Детали первопринципных расчетов . . . . . . . . . . . . . . . 46 51 52 2.2 Методика построение потенциалов межатомных взаимодействий . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.3 Техника безопасности . . . . . . . . . . . . . . . . . . . . . . 58 2

3 Результаты исследования и их обсуждение 3.1 Тестирование PAW потенциалов . . . . . . . . . . . . . . . . 60 60 3.2 Подгоночная база данных, для проверки на сходимость характеристик структур углерода . . . . . . . . . . . . . . . . . 65 3.3 Расчет сходимости в описании свойств алмаза и графита . . 66 3.4 Расширенная подгоночная база данных, для потенциала . . 75 3.5 Тестирование межатомного потенциала . . . . . . . . . . . . 77 Заключение 82 Литература 83 3

Введение Углерод является одним из наиболее широко искользуемых элементов в ядерной и термоядерной энергетике, не только в однокомпонентных системах, но и как легирующий элемент в сплавах металлических систем. Взаимное расположение атомов углерода является важным фактором, определяющим какими ствойствами будет обладать материал на его основе. Также распределение атомов углерода в твердых растворах металлов, часто имеет значительное влияние на их свойства. Чтобы иметь возможность установить механизмы процессов в условиях радиационного облучения, контролировать синтез материалов и влиять на их свойства, необходимо понимать структурное поведение, связанное с эволюцией атомов, из которых они состоят. Один из наиболее эффективных способов исследования влияния поведение атомов на свойства материалов, заключается в использовании атомистического метода моделирования, основанного на межатомных потенциалах. Первые межатомные потенциалы использовались, как иструмент, для изучения радиационных эффектов[1]. Один из методов, заслуживших высокую популярность среди межатомных потенциалов, был построен Баскесом и известен, как метод встроенного атома[2]. Из значительных недостатков данного метода, являлось, то что он не учитывал ковалентый тип связи, так как при таком типе связи преобладает анизотропный характер распределение электроной плотности. Поэтому появился модифицированный метод встроенного атома, в котором ковалентная связь учитывалась с помощью трехчастичных взаимодействий[3], причем трехчастичные взаимодействия были включены в функциональную форму межатомных не явно. Примером явного влючения трехчастичных взаимодействия явялеются потенциалы построеные Стиллинджером-Вебером. Потенциал 4

Стиллинджера-Вебера получил широкое применение при моделирования жидкого кремения[4] и разрушения в графене и углеродных нанотрубках, но лежащая в основе его модели угловая зависимость, жестко связывала атом с его соседями, что делала потенциал не преносимым на другие структуры. Большой прорыв случился, когда Терсофф, построил новый потенциал, для углерода [5]. Этот потенциал учитывал ковалентностный характер связи, также включал трехчастичные взаимодействия и был переносим на большее колличество структур, благодаря тому что индивидуально атом не имел жесткой связки с однозначным набором атомов-соседей и при этом являлся зависимым от окружающей сруды в которой он находится. Вскоре данный подход был развит Бреннером, для описания углеводородов[6] и имел большую популярность. Он хорошо описывал химические реакции и учитывал кратность изменения ковалентности связей. Но ни один из этих потенциалов не включал дальнодествующие силы Ван-дер-Вальса. Что мотивировало создание новых подходов к построению межатомных потенциалов, для ковалентных систем. Сейчас развитие межатомных потенциалов продолжает являтся популярной темой в материаловедении. Очень часто это обусловленно тем, что межатомные потенциалы позволяют объяснить механизмы процессов на атомном уровне, недоступные экспериментальным методам, а также снизить стоимость эксперимента, возникающей из-за возможной высокой сложности его проведения. Большой интерес развития межатомных потенциалов ковалентных систем, возникает из-за уникальных радиационных свойств графита[7] , алмаза[8], тепловых свойств графена[9], а также исключительных радиационных свойств, материалов в составе которых есть углерод, это стали, карбиды кремния, углеродные нанокомпозиты и др. Также мотивацией для оптимизации и создания новых методов построения межатомных потенциалов, яваляется неспособность переносимости имеюимхся потенциалов на другие структуры углерода. Во многом эти проблемы возникают из-за конкурирующих процессов гибридизации и дальнодействующиъ сил взаимодействия, связанные с п-электронами. Эти проблемы пытаются решать разными способами. Так, межатомный потенциал, построенный на основе модиффицированного метода встроенного атома(ММПА)[10], даль5

нодействующие силы ван-дер-вальса учитывает, через парный потенциал Леонарда джонса. И все же он дает большую неточность в описании констатнт упругости, для графита, в виде отрицательного знака при Сij, в сравнении с положительным значением из экспериментальных исследований. Недавно был разработан новый метод построения межатомных потенциалов и опробирован на тугоплавких ОЦК-d металлах V, Nb, Mo, Ta и W [11], где ковалентная составляющая межатомной связи сильно выражена. Было показано, что точное задание трехчастичных взаимодействий через увеличение числа оптимизируемых базисных функций, описывающих взаимодействия от расстояния между атомами , позволяет воспроизводит тепловые характеристики ОЦК d-металлы, лучше чем существующие аналоги межатомных потенциалов. Целью данного исследования является определить, как влияет изменения колличества базисных фукнций в трехчастичных взаимодействиях на точность прогноза характеристик структур углерода. Для достижения цели были поставлены и решены следующие задачи: 1. Составить базу экспериментальных и полученных на основе теории функционала плотности данных; 2. В рамках нового подхода построить потенциалы для углерода, с одной, двумя, тремя и четырьмя базисными функциями в трехчастичных взаимодействиях; 3. Провести сравнительный анализ влияния изменения колличества базисных фукнций в трехчастичных взаимодействиях на точность в описнании констант упругости, для алмаза и графита. 4. Выполнить тестирование выбранного межатомного потенциала на прогноз энергии образования точечных дефектов в алмазе, энергии сублимации и диаметра молекулы фуллерена С60. 6

Глава 1 Обзор литературы 1.1 Теория функционала электронной плотности Одним из наиболее важных прорывов в физике в области компьютерного моделирование является открытие теории функционала плотности (ТФП). Теория функционала плотности представляет собой вычислительный метод используемый, для квантовомеханического моделирования в физике, химии и материаловедении, при исследования электронной структуры многочастичных систем, в частности атомов, молекул и конденсированных фаз. Популярность ТФП растет, что отражает зависимость колличества публикаций использующих ТФП в своих исследованиях представленная на рисунке 1.1 Такую тенденцию ТФП заслуживет тем, что является наиболее физически обоснованым, в сравнении другими вычислительными методами и не требует подгоночных параметров. Теория функционала плотности имеет свои корни из модели ТомасаФерми для электронной структуры материалов, ДПФ была впервые поставлена на теоретическую основу Уолтера Коном и Пьером Хоэнбергом в рамках двух теорем Хоэнберга-Кона. ДФП имеет ряд приближений, в рамках которых возможно решение уравнений Шреденгера . Уравнение Шрёдингера является постулатом, поулченное в результате обобщения экспериментальных данных[1]. Оно позволяет описать квантовой характер взаимодействующих атомов, с помощью ввода волновых функций описывающие 7

Рис. 1.1: Публикации ТФП 1985-2015 года составляющие атома: ядра и электроны. Решение уравнений Шреденгера, могло бы позволить получить практически всю всю фундаментальную информацию о свойствах материалов. Вид уравнения Шредингера представлено в формуле 1.1 ĤΨ = EΨ, (1.1) где Ĥ - гамильтаниан, представляющий полную энергию системы, Ψ и E - обозначают собственные функции и значения оператора. Но при попытки описать многочастичную систему возникает проблема сязаная с необходимостью использования большого числа волновых функций (имеющую больше двух атомов). Эта проблема в физике носит название "Многочастичная проблема"и на сегодняшний день вычисление значение гамильтониана невозможно. Вид опретора (гамильтониана), для системы взаимо- 8

действующих ядер и электронов представлен в формуле 1.2 }2 X 2 X e2 ZI e2 X 1 }2 X ∇2 e2 X ZI ZJ Ĥ = − + − + . ∇+ 2Me i |ri − RI | 2 |ri − rj | 2 MI 2 |RI − RJ | I,i i6=j I I6=J (1.2) Гамильтониан состоит из пяти операторов. Первый предсталяет кинетическую энергию электронов, второй - взаимодействие электронов и ядер, вызваное кулоновскими силами, третий - взаимодействие между электронами, вызваное кулоновскими силами, четвертый - энергия движения (кинетическая) ядер, и, наконец, пятый - взаимодействие между ядрами, вызваное кулоновскими силами. Координаты ядер обозначают прописные индексы, строчными индексами отмечены координаты электронов. MI и ZI - масса и заряд ядер I , e - заряд электрона, Me - масса электрона, r и R - положения электронов и ядер, соответственно. И все же квантово-механические взаимодействия исследуются, работы ведутся, публикуются статьи. Это возможно благодаря тому, что методы в таких работах имеют ряд приближений, которые включает теория функционала плотности. Одним из наиболее популярных является приближение, (1.2), является адиабатическое приближение. Оно значительно упрощает гамильтониан. Суть его состоит в том, что ядра движутся из-за своей массы, на несколько порядков медленнее, чем электроны, порядка скорости звука (∼ 103 м/с) , в связи с чем можно отдельно рассматривать движение электронов и ядер. Потенциал электрона по отношению к ядрам обозначают V̂ext , для других электронов V̂ . Адиабатическое приближение имеет и другое название приближение Борна — Оппенгеймера. Таким образом гамильтониан (1.2) сводится к следующему виду: Ĥ = T̂ + V̂ + V̂ext , (1.3) где T̂ - оператор кинетической энергии. Существует множество сложных методов решения многочастичного уравнения Шредингера, основанного на расширении волновой функции в детерминантах Слейтера. В то время как самым простым является метод Хартри-Фока, более сложные подходы обычно классифицируются как методы пост-Хартри-Фока. Однако проблема с этими методами заключается в огромных вычислительных усилиях, 9

что делает практически невозможным их эффективное применение в более крупных и сложных системах. Приближение Борна-Оппенгеймера частично помогает упростить решение уравнения Шреденгера, но все еще остается проблема для сложных систем имеющих колличество атомов, больших двух. Так как волновая функция зависит как минимум от 3N переменных для системы из N взаимодействующих электронов. Тем не менее, несмотря на приближение Борна — Оппенгеймера , решение уравнения все еще остается крайне сложным для систем, состоящих из более чем нескольких атомов, многочастичная волновая функция зависит как минимум от 3N переменных для системы из N взаимодействующих электронов. Помогает решить данную проблему теоремы Хоэнберга и Кона. Основная мысль превой теоремы заключается в том, что во внешнем потенциале V̂ext , для произвольной системы взаимодействующих частиц, потенциал V̂ext устанавливается однозначно. Гамильтониан (1.3) зависит от потенциала V̂ext , из чего можно заключить, что все свойства многочастичной системы можно однозначно определить, зная ρ(r). Суть второй теоремы заключается в том, что для энергии E[ρ(r)] многоцелевой функционал может быть определен в терминах ρ(r), который может быть применим для любого внешнего потенциала V̂ext . При этом для произвольно выбранного потенциала V̂ext точное значение энергии базового состояния преставляет глобальный минимум этого функционала, и плотность ρ(r), минимизирующая функционал равна точной плотности базового состояния ρ0 (r). Отсюда следует, что вторая теорема является вариационным принципом для функционала плотности E[ρ0 (r)] 6 E[ρ(r)]. Математическое описание вариационного принципа представлено в уравнении 1.4. E[ρ0 (r)] = hψ(r)|T̂ + V̂ |ψ(r)i + hψ(r)|V̂ext |ψ(r)i Z = vHK [ρ(r)] + ρ(r)Uext (r)dr, (1.4) (1.5) где ψ(r) - электронная волновая функция, Uext - внешний потенциал vHK [ρ(r)] - функционал Хоэнберга-Кона. Функционал Хоэнберга-Кона не яляется решенным, для описания квантово-механических взаимодействий. Решение этой проблемы было предложено Коном и Шэмом, они ввели дополнитель10

ную энергию, которая устанавливается разницей между энергиями реальной системы и системы независимых электронов, известаня как обменнокорреляционная энергия Exc , Функицонал Хоэнберга-Кона определяется соотношением представленным в уравнении 1.6. vHK [ρ(r)] = Ekin [ρ(r)] + EHF [ρ(r)] + Exc [ρ(r)], (1.6) где Ekin [ρ(r)] - энергия движения, EHF [ρ(r)] - кулоновская энергия взаимодействия Хартри-Фока. Однозначного вида обменно-корреляционной энергии Exc [ρ(r)] на данный момент не существует, но есть приближения в рамках которых вид обменно-корреляционной энергии может быть определен. Во многих случая именно с этими приближениями связывают ошибки электронной плотности. История развития теории функицонала плотности относительно точности определения электронной плотности представлена на рисунке 1.2 . Благодаря этому приближению теория функционала плот- Рис. 1.2: История развития теории функицонала плотности относительно точности определения электронной плотности ности же предлагает привлекательную альтернативу, будучи гораздо более универсальной, поскольку она дает возможность систематически отображать проблему многих тел для V̂ext заменяет на проблему с одним телом 11

без V̂ext . Это финальное приближение позволяет сделать расчеты в рамках теории функционала плотности осуществимыми. Одночастичный гамильтаниан предствлен в уравнении 1.7 Ĥ1 = T̂∅ + V̂HT + V̂xc + V̂ext , (1.7) где T̂∅ - оператор кинетической энергии движения электронов в изолированном состоянии. Собственные значения гамильтониана для системы из N частиц, определяются так же как и для одночастичной системы, взятой в наиболе низких φi состояних, для всех N. Поскольку частицы в системе Кона-Шама являются невзаимодействующими фермионами, волновая функция Кона-Шама является единственным детерминантом Слэтера, построенным из набора орбиталей, которые являются наименьшими энергетическими решениями этих уравнений. Уравнение Кона-Шема представлены в 1.8 Ĥ1 φi (r) = Ei φi (r). (1.8) Эти уравнения являются базовыми для вычислительного метода теории функционала плотности. Уравнения Кона-Шэма [12] решаются путем разложения одноэлектронных волновых функций по базису и диагонализации матрицы гамильтониана. Полученные орбитали используются для расчета электронной плотности основного состояния системы взаимодействующих частиц: ρ(r) = N X φ∗i (r) × φi (r). (1.9) i=1 Теоремы Хоэнберга и Кона гарантируют, что электронная плотность основного состояния однозначно определяет гамильтониан, поэтому гамильтониан Кона-Шэма выражается через плотность ρ. Сумма таких орбитальных энергий является полной энергией ситемы. Z N X δExc [ρ] ρ(r)dr. E= i − EH [ρ] − δρ(r) i=1 (1.10) Задача решается при нахождении начальной зарядовой плотностии , путем подстановки в уравнение выбраного значения и самосогласованных решений уравнения Кона-Шэма путем его итераций. Основным камнем предкновения при решении уравнения Кона-Шема является обменно-корреляционный 12

функционал, от него зависит не только точность, но и вычислительная стоимость расчетов ТФП. Функционал может быть выражен следующим образом 1.11: Z Êxc [ρ↑ (r), ρ↓ (r)] = ρ(r)xc [ρ↑ (r), ρ↓ (r)]dr, (1.11) где подинтегральная функция представляет плотность энергии, и xc [ρ↑ (r), ρ↓ (r)] является обменно-корреляционной энергией приходящийся на один электрон. ρ↑ (r) и ρ↓ (r) являются электронными плотностями со спином вверх и спином вниз, соответственно. В случии, если спиновая-поляризация отсутствует , необходима только одна зарядовая плотность. Нахождения xc [ρ↑ (r), ρ↓ (r)] может быть реализованно используя приближения с различными уровнями сложности. Приближение локальной (спиновой) плотности, является одним из самых простых. В нем приближение для обменно-корреляционной энергии, зависящее только от плотности, может развиваться по-разному. Наиболее успешный подход основан на однородном электронном газе. Он строется путем размещения N взаимодействующих электронов в объеме V , с положительным зарядом ядра, поддерживающим систему нейтральной. Затем N и V переносятся на бесконечность таким образом, чтобы плотность ρ = N/V ) была конечной. Это полезное приближение, так как полная энергия состоит из вкладов только от энергии кинетической энергии и обменной корреляции и что волновая функция выражается в терминах плановых волн. В частности, для постоянной плотности ρ плотность обменной энергии пропорциональна ρ1/3 . Многие подходы могут давать локальные аппроксимации обменно-корреляционной энергии. Однако в подавляющем большинстве успешные локальные аппроксимации являются теми, которые были получены из модели гомогенного электронного газа (HEG). В связи с этим LDA обычно является синонимом функционалов, основанных на приближении HEG, которые затем применяются к реалистичным системам (молекулам и твердым телам. В рамках L(S)DA величина xc [ρ↑ (r), ρ↓ (r)] представляет собой обменно-корреляционную энергию с однородными спиновым плотностями, приходящаяся на частицу электронного газа ρ↑ (r) и ρ↓ (r). Эта энергия была точно найдена в работе [13], и имела несколько параметризаций которые и в настоящее время доступны, 13

для работы с программами, реализующими теорию функционала плотности. Значение обменно-корреляционной энергии в рамках L(S)DA является точным для спиновых плотностей, медленно меняющихся по пространству, но было показано, что оно может быть также точным для многих твердых тел с быстрыми вариациями плотности. L(S)DA является лучшим приближением для обменно-корреляционной энергии, чем обменная и корреляционная энергия по отдельности, потому что эти две энергии имеют противодействующие нелокальности, которые стремятся взаимно аннулировать друг друга. Применимо к задачам квантовой химии это приводит к корректному предсказанию геометрий и колебательных частот. Значительным недостатком этого приближения является неверный прогноз немагнитного и антиферромагнитного состояний чистого железа, как более стабильных по сравнению с ферромагнитным состоянием. Если учитывать градиенты электронных плотностей, обменно-корреляционный функционал будет иметь дополнительные переменные, xc [ρ↑ (r), ρ↓ (), ∇ρ↑ (r), ∇ρ↓ (r)]. Для спиновых плотностей, которые медленно изменяются в пространстве, главная поправка к L(S)DA приводит к приближению коррекции градиента второго порядка (2GA), но оно дает серьезную избыточную коррекцию для реалистичных плотностей. В формализме Кона-Шэма каждый электрон в системе эффективно окружен обменно-корреляционной дыркой, которая приводит к уменьшению его потенциальной энергии. Это результат сохранения числа электронов, принципа запрета Паули и кулоновского отталкивания. 2GA избыточно корректирует плотность, потому что оно нарушает некоторые правила дырок. Лангрет и Пердью обнаружили, что корреляционная дырка нарушает правило суммирования, которое гарантирует, что Кулоновское отталкивание не меняет интеграл заряда дырки даже если оно меняет конфигурацию дырки. Пердью также обнаружил, что плотность обменных дырок не везде отрицательная. Это означает, что дырка может существовать там, где электрон отсутствует. Это побочное явление случается относительно далеко от центра электрона, и обобщенное градиентное приближение [14, 15] параметризуется таким образом, чтобы соответствующие обрезания могли восстановить свойства дырок. Приближение GGA с параметризацией 14

PBE (Perdew-Burke-Erzenhorf) [], используемое в данной работе, достигает тех же результатов без введения каких-либо параметров. Обобщенные градиентные аппроксимации (GGA) выходят за пределы описаний LDA и LSD, включая градиенты плотности и значительно улучшают расчетные результаты Существуют также другие приближения, которые в некоторых случаях могут давать более точные результаты по сравнению с L(S)DA и GGA, однако их использование зачастую увеличивает время расчетов в разы. В нашей работе не было необходимости использовать такие приближения, поэтому их описание здесь приводиться не будет. 1.1.1 Псевдопотенциалы В вычислительных методых основанных на теории функционала плотости также в качестве приближения для упрощенного описания сложных систем используется псевдопотенциал или эффективный потенциал. Псевдопотенциал представляет собой попытку заменить движение электронов и ядер на движение только валентных электронов и оставшегося ионного остова. Это позволяет описывать псевдоволновые функции с гораздо меньшим числом мод Фурье, что делает практические применения плоских волн, менее вычилительно затратным. В этом подходе обычно рассматриваются только химически активные валентные электроны, а электроны ядра «заморожены», рассматриваемые вместе с ядрами как жесткие неполяризуемые ионные ядра. Можно самосогласовать псевдопотенциал с химической средой, в которую он встроен, что оказывает влияние на ослабление ионного остова, хотя это делается крайне редко. Псевдопотенциалы являются производным от атомного эталонного состояния, требуя, чтобы псевдо- и все-электроны валентного состояния имели одну и ту же энергию и амплитуду (и, следовательно, плотность) вне выбраного радиуса отсечки. Псевдопотенциалы с большим радиусом отсечения считаются более мягкими, то есть более быстрыми сходящимися, но в то же время менее переносимыми, то есть менее точными, чтобы воспроизводить реалистичные объекты в разных средах. В данной работе используются псевдопотеницалы, которые были построенны в рамках метода спроецированных плоских волн (PAW 15

- projected augmeneted wave) [16]. Метод расширенной волны проектора (PAW) является обобщением методов псевдопотенциала и линейной расширенной плоской волны, в результате чего позволяет проводить вычисления теории функционала плотности с большей вычислительной эффективностью. Валентные волновые функции имеют тенденцию к быстрым колебаниям вблизи ионных остовов из-за того, что они ортогональны состояниям ядра. Эта ситуация является проблематичной, поскольку она требует много компонентов Фурье, для точного описания волновых функций. Подход PAW решает эту проблему, трансформируя эти быстро осциллирующие волновые функции в гладкие волновые функции, которые являются более удобными для вычислений и обеспечивают способ вычисления всех электронных свойств этих гладких волновых функций. Валентную волновую функцию можно получить в результате линейного преобразования псевдоволновой функции T̂ : |ψi = T̂ |ψ̃i. (1.12) Такое преобразование возможно только только внутри сферы, считая за центр ядро атомов, так как здесь волновая функция ведет себя подобно осциллятору. Если рассматривать валентную волновую функцию и псевдоволновую за пределами за пределами сферы, то их поведение юудет индентичным. Таким образом, гладкая псевдоволновая функция находящаяся внутри сферы может быть разложена по сферическим гармоникам |φ̃i i: |ψ̃i = X |φ̃i ici . (1.13) i Предложенный тип разложения очень удобен, так как |φ̃i i могут быть взяты при решении, для свободного атома уравнения Кона-Шэма. Учитывая, что все преобразования являются линейными T̂ , коэффициенты ci будут определятся как: ci = hp̃i |ψ̃i, (1.14) где p̃i - функции проектора, которые удовлетворяют условию: hp̃i |φ̃j i = δij , 16 (1.15)

где δij = 1 при i = j, и δij = 0 при i 6= j. Волновая функция |ψi может быть разложена внутри сферы таким же способом как |φ̃i i: X |ψi = |φi ici , (1.16) i и тогда полная волновая функция может быть выражена следующим образом: |ψi = |ψ̃i + X |φi ici − i X |φ̃i ici . (1.17) i Здесь полную волновую функцию |ψi можно получить из псевдоволновой функции |ψ̃i если в области ядра заменить псевдоволновую функцию P P i |φ̃i ici точной волновой функцией i |φi ici . Таким образом, PAW метод используя разложение волновых функций в плосковолновой базис и затем применяя набор проекторов на радиальные сетки в центрах атомов для того, учитывает в расчетах волновые функции ядерных электронов. Уникальность этого метода заключается в том что он сохраняет полную волновую функцию и зарядовую плотность всех электронов и при этом остается эффективным с точки зрения ресурсных затрат на расчеты. Ключевым параметром псевдопотенциала является то, каким способом было произведено разделение валентных и остовных электронов. Для атомов это разделение можно сделать основываясь на их электронной конфигурации. Например для углерода, имеющего электронную оболочку 1s2 2s2 2p2 , в качестве остовных можно выбрать 1s2 электроны, оставив в качестве валентных только 2s2 2p2 электроны. 1.1.2 Расчеты электронной структуры В кристалографии при анализе таких явлений, как дифракция, рассеяние фотонов, движение электронов в потенциальном поле, связанных с периодичностью расположения частиц, значимую роль играет введение математической абстракции известной, как обратная решетка, позволяющая математически просто описать условия протекания того или иного явления в твердом кристаллическом теле. Уравнений Кона-Шэма решаемое методом самосогласования, может быть проведено как в прямом, так и в обратном пространстве. Однако, следует заметить, что в обратном пространстве 17

решить данную задачу во значительно проще. К примеру, одной из таких задач явялется задача расчета электронной структуры твердого тела. Рассмотрим некоторые аспекты численной реализации метода функционала плотности. Собственные состояния произвольного уравнения Шредингера, для независимой частицы, характеризующие движение электрона в эффективном потенциале V̂ef f (r), удовлетвлетворяют уравнению на собственные значения: }2 2 ∇ + V̂ef f ψi (r) = Ei ψi (r), Ĥef f (r)ψi (r) = − 2Me (1.18) где Ĥef f )(r) - эффективный гамильтониан. Уравнения Кона-Шэма (1.8) очевидно удовлетворяют приведенному выше уравнению. Введение граничных условий Борна-Кармана в уравнение (1.18) сделает волновую функцию ψi периодичной, так что она может быть разложена в полный набор компонентов Фурье: ψi (r) = X q X 1 ζi,q × √ exp(iqṙ) ≡ ζi,q × |qi. V q (1.19) Плоские волны ортонормированные, поэтому hq0 |qi = δq,q0 . В уравнении (1.19) V - объем рассчетной ячейки, транслируемой в пространстве, ζi,q коэффициент разложения волновой функции ψi по базису ортонормированных плоских волн |qi. Запишем представление уравнения Шредингера в обратной решетке: X hq0 |Ĥef f |qiζi,q = Ei X q hq0 |qiζi,q = Ei ζi,q0 . (1.20) q Потенциал V̂ef f (r) также периодичный и имеет представление в обратной решетке. Матричные элементы потенциала в обратной решетке следующие: 0 hq |V̂ef f |qiζi,q = X V̂ef f (Gh )δq0 −q,Gh . (1.21) m В этом случае символ Кронекера δq0 −q,Gh имеет ненулевые значения только тогда, когда q и q0 отличаются на некоторый вектор обратной решетки Gh . Обозначение можно сделать проще, если ввести определения q = k + Gh и q0 = k + G0h . Оператор кинетической энергии (}2 /2Me )|q|2 δq,q0 , так что 18

уравнение Шредингера для произвольно заданного k можно зафиксировать в матричной форме, представленной в уравнении1.22: X Ĥh,h0 (k)ζi,h0 (k) = Ei (k)ζi,h (k). (1.22) h0 По симметрии k-точки внутри зоны Бриллюэна являются неэквивалентными, поэтому только эти точки требуют дополнительного рассмотрения. Важно отметить, что плотность k-точек в обратном пространстве с увеличением объема расчетной ячейки растет и стремится к непрерывной плотности. Для любой k существует свой дискретный набор состояний i = 1, 2, . . . с определенными значениями Ei , которые становятся непрерывными зонами. Набор энергетических зон и k-точек представляет собой электронную зонную структуру, оперируя полученными характеристиками этой структуры, можно получить много информации о твердом теле. Наиболее часто используемая на практике величина - плотность электронных состояний (DOS - Density of states), которую можно получить суммируя все зоны и k-точки: g(E) = 1 X δ(Ei (k) − E), Nk (1.23) i,k где Nk , полное число k-точек, является нормирующим множителем. Определив плотность электронных состояний можно вычислить электронную энтропию: Z Selec (T ) = −3kB inf g(E)[(1 − fT (E))ln(1 − fT (E)) + fT (E)lnfT (E)]dE, 0 (1.24) где fT (E) - функция распределения Ферми-Дирака: fT (E) = 1 . exp(E/kB T ) + 1 (1.25) С точки зрения вычислительных ресурсов в расчетах, наиболее тяжелой задачей, проводимой в рамках теории функционала плотности, является диагонализация матрицы гамильтониана (1.22), которая слишком огромной для реальных систем, из-за чего не может быть решена напрямую. Однако существует общий подход для решения этой задачи, он состоит 19

в итеративной диагонализации матрицы до тех пор, пока не достигнется сходимость. Сходимость определяется путем изменение параметра между смежными диагонализациями, который в свою очередь должен быть меньше, определенного значение. В качестве этого параметра может быть силы, полная энергия, собственные значения и другие. Существует множество алгоритмов, для реализации данного подхода, такие как метод сопряженного градиента, метод остаточной минимизации и другие. Также важным факторами, влияющим на производительность вычислений, являются число k-точек в зоне Бриллюэна и число плоских волн. Большинство вычислений выполняются, учитывая только k-точки в непроводимой части зоны Бриллюэна, которая является уменьшенной за счет симметрии кристаллической решетки. Например, зона Бриллюэна в кубических решетках может быть уменьшена до 1/48 от начальной. Одним из наиболее общих способов распределения k-точек в зоне Бриллюэна является способ основанный на построение сетки, представленный в уравнении 1.26: kn1 ,n2 ,n3 ≡ 3 X 2ni − Ni − 1 2Ni i Gi , (1.26) где ni = 1, 2, . . . , Ni . Представленная выше схема была предложена Монкхорстом и Паком переносима на любые структуры, используя единообразный набор k-точек. Число плоских волн, по базису которых раскладывается волновая функция ψi , определяется энергий обрезания Ecut следующим образом: }2 |G + k|2 < Ecut . (1.27) 2Me Таким образом можно заключить, что число плоских волн различается в каждой k-точке. Это объясняется тем, что число плоских волн плавно меняется с изменением объема. Такой подход значительно улучшает точность расчетов энергии и объема, которые используются, например, для нахождения равновесных параметров решетки системы. В общем случае, число k-точек и значение энергии обрезания Ecut являются одними из основных параметров, по которым проверяется сходимость полной энергии системы. Есть также и другие параметры, по которым сходимость рассчитываемых величин проверяется , например, по ширине размытия уровня 20

Ферми, размеру расчетной ячейки и другие. Целесообразней подобирать такие значения параметров сходимости, которые бы могли позволить достичь требуемой точности расчетов с минимальными затратами вычислительных ресурсов. 1.2 Теория динамики решетки Описание термодинамических свойств кристалов, можно связывать с тепловыми колебаниями атомов, которые в свою очередь описаны в рамках теории динамики решетки. Одной из популярных модейлей динамики решетки, ялвяется модель предложенная Борном и вон Карманом [] в 1912 году, которая в современном материаловедении и физики твердого тела, является одним из краеугольных камней. Данная модель представляет механический подходом к решению задачи о стабильности решетки, который в свою очередь является мощным механизмом, при изучении частот колебаний атомов решетки. В литературе существует много статей и книг, посвященных теории динамики решетки, в этой работе мы ограничемся только основными аспектами этой теории. 1.2.1 Модель Борна-вон Кармана Если предположить, что, в результате термических флуктуаций, атом k расположенный в ячейке кристалла l (атом lk) смещаетится из положения равновесия xlk за определенное время t, то такое смещение можно описать вектором ulk (t). Таким образом, можно определить положение атомов в определенный момент времени t искходя и ниже представленного уравнения 1.28: rlk (t) = xlk + ulk (t). (1.28) При адиабатическом приближении Борна-Оппенгеймера можно говорить, что в произвольный момент времени t потенциальная энергия U кристалла, должна определяется толко функцией зависящей от мгновенных положений атомов в системе. Таким образом, если расложить в ряд Тейлора потенциальную энергию кристалла в декартовых координатах по степеням 21

атомных смещений α = x, y, z, получим следующее уравнение 1.29: U = U0 + XX lk α φlk,α ulk,α + 1 XXXX φlk,l0 k0 ,αα0 ulk,α ul0 k0 ,α0 +. . . , (1.29) 2 0 0 0 α lk α lk Важно отметить, что благодаря симметрии кристалла, большое колличество смещений будут эквивалентными, что упрощает решение задачи. В предсавленном выше выражении (1.29), чтобы найти коэффициенты разложения в ряде Тейлора, необходимо взять производные U по смещениям, представленные в уравнении 1.30: φlk,α и φlk,l0 k0 ,αα0

∂U

= , ∂ulk,α

∂ 2U

= ∂ulk,α ∂ul0 k0 ,α0

0 (1.30) (1.31) Вычисления производных необхоимо проводить в точках, которые соотвествуют равновесным положениям атомов, в этом случае все смещения равны нулю, из-за того что в равновесном положении на атомы не действуют силы, соотвественно и производная (1.30) в тих точках тоже будет равна нулю. Эти производные вычисляются в точках, соответствующим равновесным положениям атомов, когда все смещения равны нулю, Производная (1.31) можно также выразить используя матричную форму, представленную в виде матрицы межатомных силовых постоянных измеренных между парами атомов lk и l0 k 0 , матрица силовых постоянных представлена в уравнении 1.32: Φlk,l0 k0   φ φ φ  xx xy xz   = φyx φyy φyz  , φzx φzy φzz (1.32) Важно отметить, что элементы матрицы Φlk,l0 k0 все также относятся к паре атомов lk и l0 k 0 , но для удобства обозначения их индексы здесь опущены. Здесь, компонента матрицы Φlk,l0 k0 - это есть сила, которая действует на атом lk в направлении α и возникает при смещении атома l0 k 0 в направлении α0 . Межатомные силовые постоянные являются важным парметром, при расчете термодинамических свойств, часто они могут быть использованы для получения уравнений движения атомов, которые в свою очередь 22

совершают могут быть представлены в виде гармонических колебаний в кристалле, уравнений движения атомов представлено ниже 1.33: X Φlk,l0 k0 ul0 ,k0 Mk ülk = − (1.33) l0 k 0 Решение уравнения (1.33) так же может быть найдено в виде плоских волн, для этого досточно выполнить периодически граничные условия ul,k = √ 1 ejk (q)exp(i[qẋl − ωj (q)t]), Mk (1.34) где q - волновой вектор, ωj (q) - круговая частота, ejk (q) - вектор поляризации, j - индекс ветви. Можно заметить, что в результате подстановке смещений, которые вызванны распространяющейся волной, в уравнения движения (1.33), ей эквивалентным становится Фурье-преобразование правой части этих уравнений. ТВ результате чего данная задача сводится к решению предстваленного ниже уравнения 1.35: X ωj2 (q)ejk (q) = Dkk0 (q)ejk0 (q), (1.35) k0 Соответствующие собственные вектора могут быть посчитаны для всех атомов k базиса решетки, эти вектора являются векторами поляризации колебаний решетки ejk0 (q). Так же можно рассчитать путем диагонализации матрицы Dkk0 (q), соответствующие им круговые частоты ωj (q). Таким образом, базис решетки состоящий из N атомов k, то D(q) будет матрицей 3N ×3N , которая построенна из 3×3 на основе подматриц Dkk0 (q), которые в свою очередь представлены Фурье-образами Φlk,l0 k0 . Для 3N собственных векторов и собственных значений динамической матрицы D(q), можно найти, что , если они вычисленные при конкретном волновом векторе q, то должны соответствовать 3N собственным модам колебаний кристалла для представленного волнового вектора. ωj (q) называются ветвями дисперсии фононов j. Каждое значение q отображает конкретную моду. Динамическая матрица D(q) является Эрмитовой для любого значения q, поэтому ее собственные значения ωj2 (q) являются действительными, а ωj (q) могут быть как полностью действительными, так и полностью мнимыми. Также значения ωj (q) могут быть мнимыми, тогда система считается механически нестабильной. 23

Заметим, что только определенные значения q моогут быть разрешены. Они определяются из условия, если они не затухают в решетке. Для того чтобы определить эти значения в кристалле, который является бесконечным, используют периодические граничные условия. Такие условия обчыно называют граничными условиями Борна-вон Кармана. Принцип граничных условий Борна-вон Кармана, заключается в следующем: бесконечный кристалл разделяются на бесконечное число конечных ячеек, в каждой из которых содержится очень большое число N примитивных ячеек. Для примера, если считать каждую конечную ячейку кубической имеющей размеры La1 , La3 , La3 , так что N = L3 . Будет выполнятся слудующее условие: ulk = ul0 k0 , (1.36) где l = l1 , l2 , l3 l0 = l10 + L, l0 = l20 + L, l0 = l30 + L. Можно заметить, что этом случае смещения атомов, будет приводить к смещению эквивалентных атомов в разных конечных ячейках, при этом они будут отличаться только на фазовый множитель. Такое предложение принято для математического удобства, иногда из-за него возникают ошибки, например, когда ячейка слишком мала, а смещения большие, но тем не менее оно является неотъемлемым в современной науку при изучении материалов, вычислительными методами на основе компьютерного моделирования. Решетка может поддерживать волны, заданные выражением (1.34), поэтому: q = η1 b1 + η2 b2 + η3 b3 , (1.37) где ηi = ni /L для i = 1, 2, 3 и ni - целые числа. Плотность q-точек будет ˙ 2 × a3 )|, объем примитивной ячейки в прямом проN v/8π 3 , где v = |a1 (a странстве. Заметим, что важной величиной при определении кривых плотности и дисперсии фононов является плотность q-точек, поскольку именно она определяет из разрешение. Плотность q-точек обычно задается также, как и плотность k-точек. 1.2.2 Метод сверхячеек При расчете частот фононов используют метод ячеек, который позволяет учитывать влияние на частоты фононов условий периодичных границ. 24

Так например, смещение одного атома вызывает аналогичное смещение всех его периодических образов. Поэтому силы, которые возникают в результате его смещения из положения равновесия, будут убывать если растояние между атомами будет расти. Что может вызвать ошибки в расчетах частот фононов. Так как частоты фононов расчитываются путем построения динамической матрицы D(q), после диагонализации. А динамическая матрица строится путем Фурье-преобразование матрицы силовых постоянных Φlk,l0 k0 . Матрицца силовых постоянны строится, используя прямой расчет силы в сверхячейке, действующие на атомы, которые возникают при смещении симметрично неэквивалентных атомов. Одним из способов расчета этих сил основан на теореме Хелльмана-Фейнмана [Feynman1939], в готорой говорится, что сила, действующая на произвольный атом в системе электронов и ядер может быть посчитана, если известны только координаты ядер атомов, взятые относительно других ядер и распределение электронной плотности. Для избежания ошибок в результате влияния атома на свои образы, строится сверхячейка, размеры которой подбираются таким образом, чтобы силы, действующие на другие атомы становились пренебрежимо малы при увеличении расстояния от образа атома до действиельного смещения. Так же важно отметить, что величина смещения не должна быть больше некторой величины, которая могла бы гарантировать сохрание гармонического основания при описании теловых колебаний. Обычно величину смещения задают равной 0.01-0.04 Å. Это связано с тем что структура кристала при рассчтах в рамках теории функционала плотности релаксируется до нибольшей величины сил, действующей между атомами. Поэтому смещение межу атомами, должно вызывать силы на один или три порядка большие, чем в полностью отрелаксированной структуре. При расчете частот фононов строят динамическую матрицу D(q), затем проводят для нее диагонализацию. Модель сверхячейки для расчета силовых постоянных требуется потому, что при расчетах в рамках теории функционала плотности применяются периодические граничные условия и смещение одного атома будет означать смещение всех его периодических образов. 25

1.2.3 Вклад тепловых колебаний в термодинамические свойства твердых тел Термодинамические свойства системы моут быть полностью определены частотой фононов ωj (q), отвечающие за теловые колебания атомов. Так, например, можно найти внутреннию энергию системы атомов, совершающих гармонические колебания, используя уравнение 1.38, представленное ниже: E= 3N N X X q }ωj (q)[nj q + 1/2], (1.38) j=1 где nj q = 1/(exp(}ωj (q)/kT ) − 1) - число фононов с частотой ωj (q). Так, если фононы отсутствуют, то их энергия (1.38) носит названиие энергии нулевых колебаний. Свободная энергия Гельмгольца есть: F = −kT lnZ, (1.39) где Z - статистическая сумма, которая выражается следующим образом: }ωj (q) Y exp − 2kT . (1.40) Z= }ωj (q) qj 1 − exp − kT Подставляя (1.40) в выражение для свободной энергии Гельмгольца получим: F = kT X }ωj (q) qj 2kT + ln[1 − exp(−}ωj (q)/kT )]. (1.41) Через свободную энергию Гельмгольца можно получить энтропию тепловых колебаний, использую стандартное термодинамическое соотношение: ∂F (1.42) ∂T X 1 X = }ω(q, j)coth(}ω(q, j)/2kT )−k ln[2sinh(}ω(q, j)/2kT )] (1.43) 2T q,j q,j S=− Термодинамические функции выраженные через частоты фононов, имеют множество приложений. Здесь мы ограничемся, приведя только три основных их применения: Первое, они помогают понять механизмы процессов термадинаических вкладов колебаний атомов, на уровне детализации 26

обосновать ермодинамические свойства твердых тел полученные экспериментально. Второе, полученные результаты могут представлять болшую ценность в качестве прогнозов термодинамических характеристик материалов из-за высокой стоймости эксперимента или невозможности его проведения. И, наконец, эти результаты могут являются уникальной находкой, как приложение, позволяющее объяснить стабильности кристаллов и фазовых переходов. 1.3 Молекулярная динамика Получить информацию о динамической эволюции крупномоштабных систем атомов рассматриваемой в течении фиксированного периода времени стало возможным благодаря развитию молекулярной динамики, которая является вычислительным методом для изучения динамики систем взаимодействующих атомов и во многих случаях основана на межатомных потенциалах [17]. В модели молекулярной динамики движения атомов имтируется на основе уравнений движений Ньютона, как функция от времени. С помощью межатомных потенциалов вычисляется потенциальная энергия и силы взаимодействующих атомов, после, методом молекулярной динаимки, определяются траектории движения атомов, путем численного решения уравнений Ньютона. Согласно закону Ньютона сила действующая на атом прямопропорциальна его ускорению и имеет вид представленый в выражении (1). dr ri (t) = Fi (t), (i = 1, 2.., N ) (1.44) mi dt2 где Fi - сила действующая на атом i, m и ri - масса и положение атома i, соответственно. Fi = −∇i V (r1 , r2 , ..., rN ) (1.45) где V (r1 , r2 , .., rN ) потенциальная энергия системы, которая зависит от положения атомов в системе. Для сложения уравнений движения Ньютона применялся численный метод - алгоритм Верле в скоростной форме [18]. Данный метод часто используется для расчёта траекторий частиц в молекулярно-динамическом моделировании. В предложенном методе ско27

рости и положения частиц вычисляются при одном и том же заданном значении времени. Поэтапное работа алгоритма Верле [18] представлена в скоростной форме и имеет следующий вид : 1. Значение скоростей частиц можно получить из следующего уравнения: ~v (t + 12 ∆t) = ~v (t) + 12~a(t)∆t, 2. Координаты частиц определяют следующим уравнением: ~x(t + ∆t) = ~x(t) + ~v (t + 21 ∆t)∆t, 3. Для расчёта ускорения частиц уравнение ~a(t + ∆t) с использованием потенциалов межатомных взаимодействий при положениях частиц ~x(t + ∆t), 4. Скорости же частиц можно рассчитать используя уравнение~v (t + ∆t) = ~v (t + 12 ∆t) + 12~a(t + ∆t)∆t. Так же следует отметить, описанный выше алгоритм может быть сокращен, если исключить этап расчета скоростей на половине шага ∆t молекулярнодинамического моделирования. 1. Тогда координаты частиц определяются следующим образом: ~x(t + ∆t) = ~x(t) + ~v (t)∆t + 12~a(t)∆t2 , 2. Ускорений частиц ~a(t + ∆t) определяют с использованием межатомных потенциалов определяющих взаимодействия частиц в определных положениях ~x(t + ∆t), 3. Скорости же частиц можно рассчитать используя уравнение ~v (t + ∆t) = ~v (t) + 21 (~a + ~a(t + ∆t))∆t. При этом, согласно представленному ваше алгоритму, ускорения частиц ~a(t + ∆t) зависят только от их положений ~x(t + ∆t) и не зависят от их скоростей ~v (t + ∆t). Однако, нужно отметить, что также как и в ТФП, при построении межатомных потенциалов, также используют ряд приближений. Первое - приближение Борна-Оппенгеймера, в котором говорится, 28

что динамика электронов происходит настолько быстро, что мгновенно реагируют на движение ядер, адиабатически следуя за ними. Как следствие, их можно рассматривать отдельно. Второе относится к ядрам, которые являются более тяжелыми, чем электроны и описываются классической ньютоновской динамикой. Одним из основных достоист межатомных потенциалов является время расчетов. Так при молекулярном моделировании с использованием межатомных потенциалов, время моделирования возрастает линейно с увеличнием числа атомов. Используя ТФП,время моделирования возрастает прямо пропорционально увеличнием числа атомов в кубе. Например, если даже использовать известные программные пакеты (Abinit, VASP), которые позволяют реализовать теорию функционала плотности (ТФП) подчеркивая все современные возможности суперкомпьютеров, представляется возможным производить моделирование ячейки состоящей всего лишь из нескольких тысяч атомов. Но для изучения многих явлений таких как пластическая деформация, фазовые, рост зерен и диффузия, необходимо наиболее длительные методы моделирования, с образцами в которые для метода ТФП, в рамках ее нынешних возможностей невозможно получить. Проблемы, которые не могут быть решены в рамках метода ТФП, с успехом решает использование метода молекулярной динамики. Использование межатомных потенциалов, позволяет производить моделирование от миллиона до миллиарда атомов, и наблюдать за эволюцией всей системы в течении наносекунд. Использование межатомных потенциалов позволяет изучать на атомном уровне такие прцессы, как рост зерена, диффузия, пластичекую деформацию и т.п. Межатомные потенциалы имеют множество разновидностей, фундаментом которых, является функциональная форма, описывающая межатомные взаимодействия. Истинные же межатомные взаимодействия являются квантово-механическими по своей природе, и нет никакого известного способа, исходя из которых истинные взаимодействия, описанные уравнением Шредингера или уравнением Дирака, для электронов и ядер, могли бы быть приведены в аналитическую функциональную форму. Следовательно, надежность межатомных потенцилов зависит от их функциональ29

Рис. 1.3: This frog was uploaded via the project menu. ной формы. Одними из самых простых межатомных потенцаилов считаются парные потенциалы. Простейшей широко используемой моделью межатомного взаимодействия является потенциал Леннарда-Джонса . Потенциал Леонарда Джонса[19] σ σ ULD (r) = 4[( )12 − ( )6 ] r r (1.46) здесь - энергия сублимации, σ - расстояние, на котором межмолекулярный потенциал между двумя частицами равен нулю, r - расстояние между двумя частицами (измеренное от центра одной частицы до центра другой частицы). Следует отметить, что потенциал Леоннарда-Джонса является полуэмпирическим потенциалом, в котором слагаемое σ/r6 является точным выражением двухчастичных диполь–дипольных или дисперсионных взаимодействий, а слагаемое (σ/r)12 выбрано на основе качественных физических представлений. Потенциал Букиннгема[20] — функция потенциальной энергии, предложенная американским физиком Ричардом Букингемом, как приближение для энергии ван-дер-ваальсовых сил, альтернативное потенциалу Леннард-Джонса[21]. UB (r) = A exp(−Br) − C r6 (1.47) Два члена в правой части представляют собой отталкивание и притяжение, потому что их первые производные относительно displaystyle r r являются 30

отрицательными и положительными соответственно. UM (r) = De (e−2a(r − re ) − 2e−a(r−re ) ) (1.48) D — энергия связи, a — длина связи, a — параметр, характеризующий ширину потенциальной ямы. r0 — равновесное расстояние. У парных потенциалов есть некоторые неотъемлемые ограничения, такие как неспособность описать все три упругие константы кубических металлов. Потенциал Морзе, названный в честь физика Филиппа Морса, является удобной межатомной моделью взаимодействия потенциальной энергии двухатомной молекулы. Это лучшее приближение для колебательной структуры молекулы, чем QHO (квантовый гармонический осциллятор), поскольку оно явно включает эффекты разрыва связи, такие как существование несвязанных состояний. Это также объясняет ангармоничность реальных связей и ненулевую вероятность перехода для обертонных и комбинированных диапазонов. Потенциал Морзе также можно использовать для моделирования других взаимодействий, таких как взаимодействие между атомом и поверхностью. Благодаря своей простоте (всего три параметра фитинга), он не используется в современной спектроскопии. Однако его математическая форма вдохновила потенциал MLR (Морзе / Дальний радиус действия), который является самой популярной потенциальной энергетической функцией, используемой для подбора спектроскопических данных. UB (r) = A exp(−Br) − C r6 (1.49) В потенциалах многих тел потенциальная энергия включает в себя эффекты трех или более частиц, взаимодействующих друг с другом. [30] При симуляции с парными потенциалами существуют глобальные взаимодействия в системе, но они встречаются только через попарные члены. В потенциалах многих тел потенциальная энергия не может быть найдена суммой по парам атомов, так как эти взаимодействия вычисляются явно как комбинация членов более высокого порядка.В статистическом представлении зависимость между переменными вообще не может быть выражена с использованием только попарных произведений степеней свободы. Например, потенциал Терсоша [22], который первоначально использовался 31

для моделирования углерода, кремния и германия и с тех пор использовался для широкого круга других материалов, включает сумму по группам из трех атомов с углами между атомы являются важным фактором потенциала. Другими примерами являются метод внедренного атома [23], EDIP [24] и потенциалы приближения второго момента (TBSMA) с плотным связыванием, где электронная плотность состояний в области атома равна рассчитанная по сумме вкладов от окружающих атомов, а вклад потенциальной энергии является функцией этой суммы. 1.3.1 Термодинамические ансамбли Следует отметить, что траектории атомов должны генерироваться в заданном ансамбле в соотвествии с теми термодинамическими условиями в которых изучается система [25], Андерсена [26], Берендсена [27] и НоусХувера [28]. Термодинамическое состояние можно охарактеризовать функциями состояния (переменными состояния), которые определяют систему. Обычно используются следующие функции состояния: объём (V), число атомов (N), давление (P), температура (T), химический потенциал (µ), энтальпия (H), энергия (E) и энтропия (S). Так в молекулярной динамики применяют следующие ансамбли: NVE, NVT и NPT. Изотермически-изобарный ансамбль NPT или (NµT). В дополнение к данному термостату необходим баростат, при которых сохраняется количество вещества, давление и температура. Этот ансамбль играет важную роль в химии, поскольку химические реакции обычно проводятся в условиях постоянного давления. В микроканоническом ансамбле NVE система изолирована от изменений молей, объема и энергии. Это соответствует адиабатическому процессу без теплообмена. В каноническом NVT ансамбле сохраняется количество вещества, объем и температура. В NVT энергия эндотермических и экзотермических процессов обменивается с термостатом. Доступны различные алгоритмы термостата для добавления и удаления энергии из границ моделирования МД более или менее реалистичным 32

образом, приближаясь к каноническому ансамблю. Одним из частых источников замешательства является значение температуры в MD. Обычно у нас есть опыт с макроскопическими температурами, которые включают огромное количество частиц. Но температура является статистической величиной. Если имеется достаточно большое количество атомов, статистическая температура может быть оценена по мгновенной температуре, которая определяется путем приравнивания кинетической энергии системы к nk2BT , где n-число степеней свободы. Популярными методами контроля температуры являются термостат Ноус-Хувера, цепи Ноус-Хувера, термостат Берендсена, терморегулятор Андерсена и динамика Ланжевена. Термостат Берендсена может ввести эффект ледяного кубика, что приводит к нефизическим переводам и поворотам моделируемой системы. Нетривиально получить каноническое ансамблевое распределение конформаций и скоростей с использованием этих алгоритмов. Как это зависит от размера системы, выбора термостата, параметров термостата, временного шага и интегратора, является предметом многих статей в этой области. 1.3.2 Молекулярная статика Для нахождения минимальной энергии выбранной атомной конфигурации используют метод молекулярной статики. Данный метод предпологает демпфирование скоростей атомов в молекулярной динамики. Стандартная молекулярная динамика действует на каждом шаге эволюции системы. После завершения каждого шага, определяются действющие на атомы силы. При действии какой либо компоненты силы в противоположном направлению конкретного атома, этот атом останавливается и затем начинает движение в направлении этой силы. Молекулярная статика широко используется при исследовании энергии образования точечных дефектов и их комплексов. Энергия образования вакансии можно определить согласно формуле: N −1 E(N ), (1.50) N где E(N − 1, V ac) - энергия ячейки, содержащей N − 1 атомов и одну Efvac = E(N − 1, V ac) − 33

вакансию, E(N ) - энергия идеальной ячейки, содержащей N атомов. Для энергии образования межузельного атома предлагается следующая формула: N +1 E(N ), (1.51) N где E(N, Int) - энергия ячейки, содержащей N атомов и один межузельный EfI = E(N, Int) − атом, E(N ) - энергия идеальной ячейки, содержащей N атомов. 1.4 История развитие межатомных потенциалов углеродных систем Хотя метод квантово-механические методы обеспечивает высокую точность при расчетах, они являются дорогостоящими вычислительно и обычно их время растета растет прямопропорционально N 3 , из-за матричной диагонализации [51]. Напротив, существует множество применений в науке о углероде с большим числом атомов и более длительных периодов времени. Для этих ситуаций важна аналитическая потенциальная функциональная форма, которая является достаточно эффективной. Длинный список углеродных потенциалов представленный в таблице 1.2 демонстрирует существенный спрос на точные функциональные формы, которые могли бы моделировать углерод с приемлемыми вычислительными затратами. Межатомный потенциал Терсоффа [5] является самым простым и быстрым из углеродных потенциалов. Основываясь на элегантной привязке к порядку связи, он изменяет порядок связи и силы взаимодействующих атомов взависимости от количества соседей. Функция отсечки, используемая для подсчета соседей, заканчивается на сравнительно небольшом расстоянии 2,1 Å, промежуточном между первыми и вторыми соседями. Форма функции отсечки не была привязана к каким-либо экспериментальным или вычислительным данным и соответствует косинусной форме между 1,8 и 2,1 Å. Выбор короткого замыкания приносит значительные преимущества, поскольку вычисление потенциальной энергии приводит к относительно небольшим количествам атомов, что заначительно влияет на вычислиетль34

35 1990 1994 1090 Переменная отсечка, для маштабирования объема 1996 2005 Терсофф[?] Терсофф[30] Терсофф[31] Норланд[32] Терсофф[33] 2008 2013 ALREBO[35] SED-REBO[36] ReaxFF[41] ReaxFF 2015 Zhou [40] 2001 2000 2007 Mrovec [39] EDIP[24] 1999 Pettifor [38] ALREBO-M [37] 2015 2002 REBO-II [6] гибкая и общая функциональная форма Дальнее π отталкивание Репараметризация углерода улучшенная σ связь Аналитическая аппроксимация связывания Замена Леонарда Джанса потенциалом Морзе Добавлена функция экранирования Леонарда Джонса Расширенная база данных, учет короткой связи 1990 Добавляет образование углеводородных радикалов Увеличенная отсечка Увеличенная отсечка,силы Ван-дер-Ваальса Параметризирован для аморфного SiC Параметризирован, дефекты C в Si Обобщение на многокомпонентные системы REBO[34] EDIP ABOP Бреннер 1989 Терсофф [29] Репараметризорован, углерод 1988 Терсофф [5] Основне отличие Терсофф Год Название Семейство Pettifor Pettifor ALREBO REBO-II REBO REBO Терсофф Терсофф1988 Терсофф1988 Терсофф1989 Терсофф1990 Терсофф 1989 Терсофф1986 Терсофф1986 Построен на Таблица 1.1: Развитие потенцилов межатомного взаимодейтвия углеродных систем

ную скорость. С другой стороны, отсутствие какого-либо взаимодействия за пределами 2.1 Å означает, что между графеновыми слоями нет притяжения и отталкивания. Это имеет серьезные последствия для прогнозируемой плотности графита. Экспериментально графит гораздо менее плотный, чем алмаз, но с потенциалом Терсоффа оба материала имеют плотность около 3,5 г / см 3. Функция отсечки также имеет последствия при описании разрыва связей, так как форма и положение этой функции контролируют барьер активации при изменении типа гибридизации. Этот произвольный выбор обрезания, скорее всего, объясняет, почему точка плавления углерода с потенциалом Терсофф составляет около 6000 К [5] по сравнению с экспериментальным значением около 4300 К. Нордлунд аль. [8] сделал два значительных улучшения для потенции Терсофф- , добавляя долгосрочный привлекательный термин и увеличивая отсечку до 2,46 Å. В 2013 году были введены фукции экранирования Функции экранирования могут быть использованы для увеличения диапазона этих потенциалов. Мы изложим общую процедуру объединения функций экранирования с потенциалами порядка связи, которые не требуют каких-либо свойств потенциала. Мы используем этот подход для модификации работы Ж. Терсоффа [Phys. Rev. B39, 5566 (1989)], P. Erhart и K. Albe’s [Phys.Rev.B71, 35211 (2005)] и Т. Кумагайет аль. [Comput. Mater.Sci.39, 457 (2007)] Потенциалы Si, C и Si-C. потенциалы зависящие от среды были улучшены важность диапазона взаимодействия E= X Ei = i 1X Vij 2 i Vij = f (rij )[A exp(−λ1 rij ] − Bij exp(−λ2 rij ) (1.52) (1.53) где E общая энергия системы, Ei энергия отдельым мест, Vij энергия взаимодействия между двумя атомами i и j, rij расстояние между ними. A, B, λ1 и λ2 положительные числа. f(rij ) дополнительная функция отсечки введенная для уменьшения диапазона потенциала. Первый член в выражении (2) отвечает за отталкивание и получается в результате ортоганализации, U (r) = De (e−2a(r−re ) ) (1.54) Первый член в (2) равен отталкивающий, и интерпретируется как из-за ортогонализации- и т. д. 2 ’Второй термин интерпретируется в виде пред36

ставляющих связь. Следовательно, BiJ неявно включает порядок привязки и должен зависеть от местных условий окружающей среды Мент. 2 ’В настоящей работе вместо введения трехмерные термины для описания сил углового угла и т. д. форма (2) строго соблюдается. Al / отклонения от простой парный потенциал приписывается зависимости от 8 и от локальной атомной среды. В частности, сила сцепления BiJ для пары ij должна быть монотонно убывающая функция числа о конкурирующих облигаций, сила конкурирующих связей и косинусов углов с конкурирующими облигации. Эти три фактора были включены в следующий простой пробный потенциал: Bij = B0 exp( zij = −zij ) b X w(rik ) [ ]n × [c + exp(−d cos θijk )]−1 w(rij ) (1.55) (1.56) k6=i,j где w(r) является "голым"потенциалом связи w(r) = fc (r) exp(−λ2 r) (1.57) Здесь zij - взвешенная мера количества облигаций, конкурирующих с облигацией ij, и b определяет, насколько быстро прочность связи падает с увеличением эффективной координации. Первый член в (3b) - это просто отношение немасштабированной связи сильные стороны связей ik и ij, поднятые до степени n. Параметр n, таким образом, определяет, насколько более близкие соседи предпочитаются более отдаленным в конкуренции для образования связей. Последний член дает зависимость от угла связи, которая берется как функция cos(k), чтобы обеспечить правильное аналитическое поведение. 8 k - угол между связями ij и ik. Простую экспоненту судили до этой функции типа Ферми, но она не обеспечивала достаточной свободы для Опишем весь диапазон углов связи точно. Заметим, что эта формулировка не является симметричной, т.е. е. , Vit j . Это просто означает, что энергия, связанная с данной связью, не разделяется поровну между двумя атомами в формальном учете (1). Однако потенциал обладает всеми физически необходимыми свойствами инвариантности. Так как выражение Абелла-Терсоффа использует идеи Полинга, он может реально описать углерод-углеродные, двойные и тройные связи длины и энергии в 37

углеводородах и в твердом графите и алмазе. Однако в промежуточных ситуациях сложения предположение о взаимодействиях ближнего соседа в сочетании с суммой по атомным узлам (см. (1)] приводит к нефизическому поведению. Например, если атом углерода с тремя ближайшими соседями связан с атомом углерода с четырьмя соседями, (1) - (3) интерполировать связь так, чтобы она была промежуточной между одной и двойной связью. Однако образование двойной связи происходит из-за перекрытия несвязанных. Поскольку атом с координацией 4 не имеет свободной орбитали, перекрытие не может происходить, и связь лучше описывается как одинарная связь плюс радикальная орбиталь. Это переплетение радикалов приводит к нефизическому описанию связывания углерода в ряде обычных ситуаций. Например, образование вакансии в алмазе приводит к образованию четырех радикалов, и поэтому это переплетение радикалов позволяет увязать энергию образования вакансий в алмазе, сохраняя при этом пригодность к графиту, если это невозможно. Первоначальный потенциал реактивного эмпирического порядка связи REBO[34] был продолжением метода Терсоффа, добавляя водород и улучшая описание радикалов. Второе поколение, известное как REBO-|| [6], было выпущено двенадцать лет спустя с более подробным описанием короткодействующего связывания. Это было достигнуто путем модификации функциональной формы с использованием улучшенной базы данных для определения параметров и добавления двухгранного термина для поряка связи отвечающий за поворот π -связей. В обоих этих потенциалах функции отсечки те же, что и у потенциала Терсоффа, за исключением того, что косинусная функция находится между 1,7 Åи 2,0 Å. Соответственно, REBO и REBO-|| также сталкиваются с проблемами связынные описнием плотности графита. Наблюдается расхождения плотности графита и алмаза , поскольку потенциал REBO-|| прогнозирует, что графит в 3 раза более плотный, чем алмаз. Широкое расширение, известное как адаптивный межмолекулярный потенциал AIREBO [37], было выпущено Стюартом в 2000 []. Хотя он был опубликован этот потенциал на два года раньше, чем REBO-||, функция AIREBO основана на REBO второго поколения. Дальнодействующие силы в нем описывается с использованием формы Леннарда-Джонса в сочетании с функцией переключения для деак38

тивации длинного ранга на коротких расстояниях. Как указано в таблице 1, было разработано значительное число вариантов REBO. Это включает в себя добавление кислорода, использование функции экранирования окружающей среды для улучшения описания образования и разрушения связи и включения уравновешивания зарядов для обозначения полярных связей. Первоначальная функциональная форма потенциала REBO представлена в уравнении 1.58 Eb = XX i [VR (rij ) − B̄ij VA (rij )] (1.58) j(>i) B̄ij = (Bij + Bji )/2 + Fij (Nit , Njt , Nijconj ) (1.59) общее число соседей для атома i, определяет связь между атомами i и j , являясь частью сопряженной системы. Переопределение радикалов теперь можно зафиксировать добавлением поправки к уравнению 1.55 для связей между парами атомов которые имеют разные координаты. Как описано ниже, нелокальные эффекты также могут быть добавлены к первому приближению для учета сопряженных и неконъюгированного слеивания. Сегодня эти две новаторские статьи имеют впечатляющий след, охватывающий почти 4000 ссылок, а баланс между ними в пользу REBO, как 2: 1, 2055/5000 Потенциал взаимодействия зависимый от окружающей среды (EDIP) для углерода [24] не связан с методами Терсофф или REBO и основан на более раннем методе EDIP для кремния. EDIP для углерода был разработан для описания роста аморфных углеродных тонких пленок, где существует конкуренция между типи гибридизации связи. Соответственно, EDIP рассматривает обрезание отлично от потенциалов Терсоффа и REBO, используя обобщенный координационный термин для захвата связей путем влючения дальнего взаимодействия отталкивания и двугранного вращения. Было обращено пристальное внимание на то, чтобы правильно предсказать энергетику создания связей и из разрыв, благодаря добавлению в базу данных квантово-мехических рассчетов по переходу из графита к алмаз. Из-за этой обработки обрезания, для потенциал EDIP относительно длинное, простирается до 3,2 Å. Помимо этой точки нет никакого взаимодействия, и, следовательно, притяжение Ван-дер-Ваальса между 39

слоями графита не описывается. Несмотря на данный недостаток, моделирование с ипользованием потенциала EDIP позволяет создавать слоистые структуры, так как во многих ситуациях движущая сила для наслаивания представляет собой продольное отталкивание между sp2-связанными слоями. Еще одним недостатком является в EDIP отсутствие энергии необходиомой, для описания изолированных sp2-связей; это происходит из-за использования атомно-центрированного порядка связи, который наследуется от родительской модели, для кремния. Потенциал реактивной силы (ReaxFF) [41] был первоначально разработан для углеводородов и с тех пор был расширен, чтобы описать большое количество атомных видов, включая кислород, нитроген, сера и многие металлы. Философия развития ReaxFF значительно отличается от большинства углеродных потенциалов. Вместо того, чтобы использовать догадки о химической связи для выбора функциональной формы, разработчики ReaxFF начали с полностью общей формы для потенциала, который был сравнен с большим набором данных. Потенциал обладает многими химическими возможностями, включая порядок связей, переноса заряда, взаимодействием Ван-дер-Ваальса. В то время как семейство ReaxFF обладало значительными успехами, такими как белки, топливные элементы и катализ, твердых углеродные фазы он описывает плохо. В качестве примера, модуль всестороннего сжатия алмаза с оригинальным потенциалом ReaxFF был примерно в 100 раз выше экспериментального значения. Для решения этой проблемы в 2015 году был обновлен набор обновленных параметров ReaxFF для твердой углеродистой фазы [42]. Этот набор параметров, известный как ReaxFFC 2013. Потенциал длинного диапазона углеродных связей LCBOP [43] концептуально похож на AIREBO, сочетая описание порядка связей с дальнодействующими взаимодействиями Ван-дер-Ваальса. Однако, в то время как AIREBO был разработан путем добавления дальнего взаимодействия к ранее существовавшему потенциалу REBO, LCBOP был разработан интегрированным образом, в котором были разработаны условия порядка связи наряду с дальним взамиодействием. Эта превосходная форма параметризации исключает проблемы с функциями переключения. LCBOP был разработан с учетом твердых углеродных фаз, и были опубликованы две основные вер40

сии. В первой версии (LCBOP-I) используется короткодействующее скоординированное обрезание 2,2 Å, аналогичное тому, которое было у Tersoff и REBO, с недостатком, заключающимся в том, что диссоциация связи происходит слишком резко из-за спада обрывов. Вторая версия (LCBOP-II) предлагает значительно улучшенное описание энергии связывания путем включения среднесрочного взаимодействия между 1,7 и 4 Å. Он был успешно применен к ряду систем , таких как жидкое состояние, тройная точка углерода и рябь в графене. Общий вид потенциальной энергии, для систем Стиленжера Веббера был раложен с использованием одночастичного, двухчастиного, трехчастичного и тд. E(1..n) = X v1 (i) + X v2 (i, j) + X v3 (i, j, k) + .. + vn (1..n) (1.60) где Е общая эергия, n-это есть число частиц. v1 одночастичные, v2 двухчастичные, v3 трехчачтичне и vn n-частичные взаимоедействия. Одночастичный потенциал v обычно описывает стенку и внешние силы, которым подвержена система. Эти для случая, рассмотренного ниже, отсутствуют, поэтому разложение (2. 1), в принципе, начинается с члена парного взаимодействия rij ) σ ri rj rk v2 (ri , rj , rk ) = εf2 ( , , ) σ σ σ v2 (rij ) = εf2 ( (1.61) (1.62) . 1.5 Межатомные потенциалы в рамках нового подхода Общая потенциалльная энергия в рамках новго подхода [11] была выведена из стандартого кластерного разложение [44, 45]. Эффективная потенциальная энергия Etot для моноатомной системы из N атомов задается выражением: Etot = N X i<j V2 (ij) + N X V3 (ijk) + i<j<k N X i<j<k<l 41 V4 (ijkl) + ..., (1.63)

где начало отсчёта для Etot было выбрано таким образом, чтобы потенциальная N атомов, разнесенных на расстояния более радиуса взаимодействия друг от друга, обращалась в ноль. Отсутствие одночастичной потенциальной функции предпологает отсутствие внешнего поля. Точность уравнения (1.63) была показана в случае суммирования используя до N частичного слагаемого и имела плодотворным подход при вычислении физических величин на малых металлических кластерах [45]. Так же можно заметить, что n-частичное представление в уравнении (1.63) используется и при описании взаимодействий между атомами в многоатомных системах как методами молекулярной динамики, так и методом Монте-Карло. Однако, важно понимать, что случае уравнение (1.63) включащая некторые слагаемые (2-, 3- , 4-называемыми частичными слагаемыми),взаимодействиями между 2, 3 или 4 атомами, являются не реальными, а представлены эффективными n-частичными взаимодействиями, которые имеют оптимизируемые параметры [45]. Существует также множество других полуэмпирических потенциалов, в которых полная энергия представляется суммой двухчастичного слагаемого и одного модельного n-частичного слагаемого, где n может достигать значений 30-50. МПА представляет известный пример таких потенциалов, в котором полная энергия описывается суммой парных взаимодействий и функции погружения, где n – число соседних атомов, дающих вклад в суперпозицию атомных плотностей [46, 47, 48, 49, 50, 51]. Здесь следует подчеркнуть, что разложение полной энергии(1.63) может потенциально описывать взаимодействия в атомной системе с необходимой точность в результате увеличения числа слагаемых и повышения точности их расчета. Представим уравнение (1.63) в при относительных положений атомов в следующих терминах: Etot = N X i<j N X N N X N X X Φ(Rij )+ G3 (cos(θjik ), Rji , Rki )+ G4 (cos(θjik ), cos(θkil ), cos i i j<k j<k<l (1.64) Здесь мы убрали ограничение для индекса i который может быть меньше следующих индексов в многочастичных слагаемых. В уравнении (1.64) Φ(Rij ) – парные q взаимодействия между атомами i и j разделенные рас~i − R ~j |2 , θjik - угол между вектором R~ji = R ~j − R ~i стоянием Rij = |R 42

~ i и т. д. Введённые функции Gn определяют nи вектором R~ki = R~k − R частичные энергии Vn , данный вид функциональной формы был впервые предложен Стилинжером и Вебером [52] на примере трёхчастичных взаимодействий V3 (ijk) = G3 (ijk) + G3 (kij) + G3 (jki) (1.65) и было далее развито Мистиоритисом с соавторами на случай четырёхчастичных взаимодействий [53]. Это соотношения можно заметить (1.64) n-частичные взаимодействия Gn являются функциями R и cos(θ). В используемом в данном диссертационном исследовании подходе [54] вводятся наборы функций одного переменного, зависящие либо от аргумента R, либо от аргумента cos(θ). При этом производится разложение функций Gn по этим наборам, как наборам базисных функций, для упрощения расчётов многочастичных слагаемых. В этом случае трёхчастичные взаимодействия взаимодействия записываются в следующем виде: E3 (N ) = n3 N X N X X i g pq (cos(θjik ))f p (Rji )f q (Rki ), (1.66) j<k p,q где g pq (cos(θjik )) - коэффициенты, зависящие от углов между связями ji и ki. Уравнение (1.66) описывает трёхчастичное слагаемое с заданной точностью при соответствующем выборе числа функций f . Следует отметить, что известные потенциалы Стилинджера-Вебера или Терсоффа аналогичны выражению (1.66) в случае одной базисной функции f и специальных выборах аналитического вида g(cos(θ)) [54]. Однако точность описания трёхчастичных взаимодействий в этих потенциалах ограничена модельным видом. Аналогичных образом записываются разложения для четырёхчастичных и последующих n-частичных слагаемых по аналогии с уравнением (1.66). Например, слагаемое для четырёхчастичных взаимодействий записывается в следующем виде: E3 (N ) = n4 N X N X X i g pqs (cos(θjik ))g pqs (cos(θkil ))g pqs (cos(θlij ))f p (Rji )f q (Rki )f q (Rli j<k<l p,q,s (1.67) 43

Таким образом, Etot представляется в виде, из которого могут быть сконструированы потенциалы различного качества посредством выбора числа n-частичных слагаемых и числа базисных функций: Etot = N X Φ(Rji ) + E3 (N ) + E4 (N ) + E5 (N ) + ... (1.68) i<j В уравнении (1.68) E3 (N ) - трёхчастичное слагаемое, E3 (N ) - четырёхчастичное слагаемое и т. д. В данном диссертационном исследовании используется только трёхчастичное слагаемое для описания зависимости Etot от углов между связями Однако, учет только 2-х частичных и 3-х частичных взаимодействий недостаточен для описания межатомных взаимодействий в металлах, для которых необходим учет большого числа n-частичных слагаемых. Для металлов сходимость разложения Etot (1.63) является медленной как показано в работе [45] из анализа относительных значений вкладов n-частичных слагаемых в энергию связи кластеров металлов. Для эффективного учёта взаимодействий высоких порядков в используемом подходе применяется ЦСП. Допустимость использования ЦСП подтверждается успехом МПА при описании металлов с ГЦК структурой [55, 47, 48]. Вид МПА выражения для Etot [46, 49] может быть получен из общего выражения для Etot (1.68) в рамках ЦСП. Согласно ЦСП трёх- и другие многочастичные слагаемые не зависят от углов между связями. Таким образом произведения функций gn (cos(θ)) в уравнениях (1.66), (1.67) и (1.68) заменяются константами Cn . Также может быть введено приближение одной функции f (R) для всех многочастичных слагаемых. Тогда для случая большого числа соседей, которое обычно рассматривается в металлах при использовании МПА (например, учёт четырёх координационных сфер в ГЦК решётке даёт окружение из 54 атомов, которые формируют вклад в суперпозицию электронной плотности на центральном атоме), уравнение 1.68 сводится к виду Etot = X i<j Φ(Rji ) + N X N X i n=3 44 n−1  Cn∗  X j6=i f (Rji ) , (1.69)

где Cn∗ = 1 n−1 Cn - постоянные. Сумма по n в уравнении (1.69) является разложением в ряд Тейлора некоторой функции F (ρi ), где ρi = X f (Rji ) (1.70) j6=i Если ρi в уравнении (1.70) интерпретировать как суперпозицию электронных плотностей f (Rji ) атомов j в месте положения атома i, уравнение (1.69) сводится к известному виду МПА [46, 49, 50]: Etot = N X Φ(Rji ) + i<j N X F (ρi ) (1.71) i Стоит отметить, что хотя уравнение (1.71) формально совпадает с видом МПА, остается принципиальное различие, поскольку, электронная плотность не может принимать отрицательных значений в МПА, тогда как физические ограничения на знак функции f (Rji ) отсутствуют [54]. При применении вышеописанных приближений для n-частичных слагаемых в уравнении (1.68) при n > 3 из уравнений (1.68), (1.66) и (1.69) получается следующее выражение для Etot [54]: С использованием ряда приближений А.Г. Липницким и В.Н. Савельевым было получено следующее выражение для Etot [54]: Etot = N X Φ(Rji ) + i<j n3 N X N X X i + N X g pq (cos(θjik ))f p (Rji )f q (Rki ) + j<k p,q F (ρ̄i ), ρ̄i = i X ρ(Rji ) (1.72) j6=i В верхнем выражении (1.72) первое слагаемое описывает парные взаимодействия Φ(Rji ) между атомами i и j. Второе слагаемое описывает трёхчастичные взаимодействия с возможностью последовательного улучшения точности посредством увеличения числа n3 базисных функций f p (Rji ), аргументом которых является расстояние Rji между атомами j и i. Функции g pq (cos(θjik )) описывают зависимость Etot от угла θjik между связями пар 45

атомов ji и ki, формирующихся вокруг атома i. Третье слагаемое в формуле (??) учитывает сумму n-частичных взаимодействий при n > 3 в рамках центрально-симметричного приближения (ЦСП). Потенциальные функции Φ(Rji ), g pq (cos(θjik )), f p (Rji ), ρ(Rji ) и F (ρ̄i ) задаются в виде кубических сплайнов (сплайнов третьего порядка). Значения в узлах сплайнов являются оптимизируемыми параметрами потенциалов. Формула для расчёта силы f~m , действующей на атом m, выводится из определения ∂E ∂E ∂E ∂E =− ; ; f~m = −∇m E = − ∂rm ∂xm ∂ym ∂zm (1.73) подстановкой в это выражение потенциальной энергии E в виде (1.72) для Etot . Давление в системе атомов рассчитывается с использованием теоремы вириала для вклада в давление, связанного с межатомными взаимодействиями, и для случая (1.72) берётся из работы [54]. 1.6 Процедура оптимизации параметров потенциалов межатомных взаимодействий Минимизация целевой функции, является одной из сложнейших задач, так как расчеты потенциала соответствующих величин из оптимизируемой базы данных, проводится в 100-мерном пространстве поиска переменных (в нашем случае в роли переменных выступают значения функций в узлах сплайнов). Минимизация целевой функции есть не что иное, что минимизация суммы среднеквадратичных отклонений величин рассчитанных. Использование нами сплайн-функций, связано с активным их использованием другими авторами в связи с возможностью использование так называемого метода сопоставления сил (force-matching method), предложенного Эрколесси и Адамсем. Этот метод стал возможен за счет включения оптимизируемых параметров при использовании сплайн-функций подгоночной базы данных из ТФП атомных сил для различных конфигураций. Для минимизации целевой функции используют известыный програмный пакет 46

POTFIT, который позволяет создавать потенциалы межатомных взаимодействий путем соотношения результатов, вычисленных этими потенциалами, с набором эталонных данных, расчитанных в рамках ТФП. Также пакет POTFIT был реализован в параллельной версии, необходимой для распараллеливания, которое производится по используемой ТФП базе данных, что позволяет строить одновременно только один. На сегодняшний день в пакет также включена возможность выбора методов парных потенциалов, МПА и ММПА при построении потенциалов межатомных взаимодействий. Оптимизация параметров потенциалов проводилась путём нахождения минимума целевой функции Z [54]: Z = Zexp + Ze + Zf , (1.74) где Zexp = i Ze = pot exp 2 exp (Ai − Ai ) , Wi 2 (Aexp i ) (1.75) pot DF T 2 ) e (Ei − ηe Ei , Wi (ηe EiDF T )2 + δe2 (1.76) X X i 2 pot ηe ~DF T ~ | X f | fi − ηr fi Zf = Wi . 2 ηe DF T 2 ~ i | + δf ηr | fi (1.77) В уравнении (1.75) Zexp - сумма квадратов относительных отклонений величин, рассчитанных с использованием построенных потенциалов, Apot i от их экспериментальных значений Aexp i , помноженных на весовые множители Wiexp . В уравнении (1.76) Ze - сумма квадратов относительных отклонений энергий, рассчитанных с использованием построенных потенциалов, Eipot от масштабированных энергий, рассчитанных в рамках ТФП ηe EiDF T , помноженных на весовые множители Wie . В уравнении (1.77) Zf - сумма квадратов относительных отклонений векторов сил, действующих на атомы, рассчитанных с использованием построенных потенциалов, f~ipot от масштабированных векторов сил, действующих на атомы, рассчитанных в рамках ТФП ηe f~DF T , помноженных на весовые множители W f . Множитеηr i i ли ηe и ηr являются оптимизируемыми параметрами потенциалов. Малые 47

числа δe и δf предотвращают образование чрезмерно больших вкладов в целевую функцию от численно малых величин. Весовые множители (Wiexp , Wie , Wif ) в уравнениях (1.75)-(1.77) нормируются следующим образом: X i Wiexp + X Wie i + X Wif = 1. (1.78) i Минимизация функции Z (1.74) проводилась с использованием специального алгоритма, который производит поиск параметров на основе смоделированного отжига [56]. Представление потенциальных функций кубическими сплайнами Если значения функции h(x) известны в конечном числе точек, то могут быть использованы различные схемы интерполяции для определения значений этой функции в промежуточных точках. Полином m − 1 степени будет точно проходить через m известных значений функции, но при этом может сильно осциллировать в промежутках между значениями функций, если m велико. Та- кая интерполяция дает плохое описание функции h(x) за исключением точек в которых она уже известна. Часто лучшей альтернативой является множе- ство кусочных функций, которые дают плавную интерполяцию между двумя известными значениями функции. При этом, на эти кусочные функции накла- дываются определенные условия их поведения в граничных точках их области определения, которые обеспечивают непрерывность функции h(x) и некото- рых ее производных в точках, где точно известны значения функции h(x). Если кусочные функции предстваляют собой кубические полиномы, то такая интерполяция носит название кубической сплайн-интерполяции. В представлении функции h(x) в виде кубических сплайнов, полином ci (x) = αi + βi x + γi x2 + δi x3 определяет поведение функции h(x) между точками xi и xi+1 , и необходимы m − 1 кусочных функций ci (x) для интерполяции между m известными значениями функции h(x). В общем 4(m−1) коэффициентов должны быть определены. Приведенные ниже условия предоставляют 4(m − 1)-2 уравнения для определения неизвестных коэффициентов, а так-же гарантируют, что кубические полиномы проходят через известные значения функции h(x), и что они непрерывны вместе со своими первыми 48

и вторыми производными:   ci (xi ) = hi (xi ), 1 ≤ i < m;        c (x ) = h(xm );   m−1 m X(ω) = ci (xi+1 ) = ci+1 (xi+1 ), 1 ≤ i ≤ m − 2;      c0i (xi+1 ) = c0i+1 (xi+1 ), 1 ≤ i ≤ m − 2;     c00 (x ) = c00 (x ), 1 ≤ i ≤ m − 2 i+1 i i+1 i+1 Граничные условия, накладываемые на первый полином c 1 (x) и на последний полином c m-1(x), предоставляют последние два уравнения необходимые для нахождения всех неизвестных коэффициентов. Обычно эти граничные условия определяются обращением в ноль вторых производных на краях интерполируемой функции:  c0 (x1 ) = 0; 1 X(ω) = c00 (x ) = 0. m−1 m Также следует отметить, что из-за разных физических соображений могут потребоваться другие варианты для граничных условий. Например, силы действующие между атомами могут определятся в результате дифференциации потенциальной энергии Etot , которые в свою очередь должны обращаться в ноль за пределами диапазона взаимодействия потенциала, что, в случае задания потенциальных функций в виде кубических сплайнов, обеспечивается в результате удовлетворения требования равенства нулю первых производных на краях области определения потенциальных функций. Кубические сплайны также могут быть использованы для определения потенциальных функций, но важно отметить, что значения этих функции не известны заранее. Их можно найти, если значения функций в узлах сплайнов являются оптимизируемыми параметрами и определяются из знания экспериментальных и ТФП характеристик атомной системы. В качетсе несомненного преимущества, является то, что кубические сплайны могут в высокой точностью воспроизводить сложные функции и тем самым предоставляют гибкие средства при описании потенциалов межатомных взаимодействий, устраняя сложную задачу получения (или угадывания) явных 49

функциональных зависимостей, подходящих для описания сложных взаимодействий. Что подтверждается их широким ииспользованием при построении межатомных потенциалов различных атомных систем, включая полупроводники [8], ОЦК металлы [58,69,74], ГЦК металлы [13, 75], ГПУ металлы [59]. 50

Глава 2 Материал и методика исследования В качестве материала исседования в данной работе были выбраны следующие структуры углерода: алмаз, графит, графен, лонсдейлит, ОЦК, ГЦК, ОЦ-8, фуллерен, α-карбин и β-карбин. Их исследование проводилось с использованием PAW потенциалов и потенциалов межатомных взаимодействий. Для оптимизации потенциалов межатомных взаимодействий были взяты в качестве экспериментальной подгоночной коснтанты упругости алмаза и графита. В качестве экспериментальной подгоночной базы данных были выбраны упругие константы пиролитического графита представленые в работе [23]. Для нахождение констант упругости использовались ультразвуковые и статические методы испытаний. Наибольшей погрешностью измерения обладют слудующие костанты упругости С11, С12 и С66 -20 ГПа. В таблице 2.1 представлены экспериментальные данные выбранные, для равновесной решётки графита, в скобках указана погрешность измерений. Следует отметить, что данные экспериментальные данные характеристики структуры графита использовались и другими исследователями [57]. Экспериментальные даные для подгонки упругих констант алмаза были выбраны из работы [58]. Для нахождение констант упругости использовался спектральный анализ основаный на Брилюеновском и Рамановском рассеяниях. Наибольшей погрешностью измерения обладет костанта упругости 51

Таблица 2.1: Экспериментальные данные по равновесной решётке графита, используемые для оптимизации параметров потенциалов. a0 параметр решётки, Cij - модули упругости Свойство Значение a0 , Å 2.46 c/a 2.72 C11 , ГПа 1160(20) C12 , ГПа 180(20) C13 , ГПа 15(5) C33 , ГПа 36.50(10) C44 , ГПа 0.40(4) C66 , ГПа 440(20) С12 -10 ГПа. В таблице 2.2 представлены экспериментальные данные выбранные, для равновесной решётки алмаза, в скобках указана погрешность измерений. Энергия сублимации, для графита была взята из справочника [59]. Параметры решетки из работы [60]. Следует отметить, что данные экспериментальные данные характеристики структуры алмаза, также использовались и другими исследователями [57], [35]. 2.1 Детали первопринципных расчетов В настоящей работе расчеты полных энергий и геометрий структур систем углерода в рамках теории функционала плотности проводились с использованием пакета программ VASP [61]. Данный пакет программ используют базис плоских волн и поддерживают формализм метода PAW. Для углерода мы использовали PAW-потенциалы, построенные в рамках обобщенного градиентного приближения (GGA) с параметризацией Perdew-Burke-Ernzerhof (PBE). В праграмном пакете VASP было представлено шесть PAW потенциалов это потенциалы C, Сd, Ch, Cs, Chnr, CGW. Для этих потенциалов были рассчитаны равновенсные параметры решетки, модули всесторонненго сжатия и свободная энергия, для графита, аламаза и графена. Для расчета характеристик структуры аламаза была использо52

Таблица 2.2: Экспериментальные данные по равновесной решётке алмаза, используемые для оптимизации параметров потенциалов. Ec - энергия сублимации, a0 - параметр решётки, Cij - модули упругости Свойство Значение Ec , эВ 7.37 a0 , Å 3.57 c/a 1.0 C11 , ГПа 1080(5) C12 , ГПа 125(10) C44 , ГПа 579(5) вана ячейка из восьми атомов, с векторами трасляции 3.57 Å, 3.57 Å, 3.57 Åи наименьшим расстоянием 1.54 Å. Для расчета характеристик структуры графена была использована ячейка из двух атомов, с векторами трасляции 2.46 Å, 2.46 Å, 20 Åи наименьшим расстоянием 1.43 Å. Для расчета характеристик структуры ОЦ-8 была использована ячейка из восьми атомов, с векторами трасляции 3.46 Å, 3.46 Å, 3.46 Å. Для расчета характеристик молекулы фуллерена была использована ячейка из шестидесяти атомов, с векторами трасляции 20 Å, 20 Å, 20 Å, при этом начальный радиус молекулы фуллерена составлял 3.35 Å. Для расчета характеристик структуры ОЦК была использована ячейка из одного атома, с векторами трасляции 2.71 Å, 2.71 Å, 2.71 Å. Для расчета характеристик структуры ГЦК была использована ячейка из одного атома, с векторами трасляции 2.17 Å, 2.17 Å, 2.17 Å. Для расчета характеристик структуры ОЦ была использована ячейка из одного атома, с векторами трасляции 2.1 Å, 2.1 Å, 2.1 Å. Для расчета характеристик структуры β-карбина была использована ячейка из двух атомов, с векторами трасляции 2.5 Å, 20 Å, 20 Å. Для расчета характеристик структуры α-карбина была использована ячейка из двух атомов, с векторами трасляции 2.5 Å, 20 Å, 20 Å. На рисунке 2.1 Представлены кристаллические структуры алмаза, графита, лонсдейлита, графена и ОЦ-8 На рисунке 2.2 Представлен вид молекулы фуллерена Для определения равновесного состояния алмаза, графена, ОЦК, ГЦК, 53

Рис. 2.1: Структуры алмаза, графита, лонсдейлита, графена и ОЦ-8 / ОЦ-8 и фуллерена, провели 100 расчетов в интервале от сжатой решетки в данную на 15 % и растянутой на 15 %. Полученные значения аппроксимировались уравнением состояния Vinet. Минимум на полученой кривой данного уравнения соотвествует, минимальному значению энергии и параметра решетки, для данного потенциала. Графит был рассчитан с равновесными параметрами решетки, взятыми из экспериментальных работ. Для определения равновесного состояния лонсдейлита, α-карбина и β-карбина, провели 100 расчетов в интервале от сжатой решетки на 10 % и растянутой на 10 %. Из полученных значений собрались изолинии. Минимум на получнных картинках определялся путем перемещения сетки, состоящей из 10000 на 10000 узлов, для поиска наиболее яркого участка изолиний, который и соотвествует минимальному значению энергии и параметра решетки, для данного потенциала. Для всех предложенных структур значения Ecut и плотности k-точек, подобранных для чистого углерода, было достаточно, чтобы сохранить заданную точность расчетов полной энергии. При этом сетка строилась таким образом, чтобы плотность k-точек во всех направлениях была примерно одинаковой. Размытие уровня Ферми было получено методом Methfessel-Paxton [62], а ширина размытия составляла 54

Рис. 2.2: Молекула фуллерана 0 эВ, при расчетах в пакете VASP. Основным критерием сходимости расчетов по параметру Ecut и наборам k-точек в программном пакете VASP являлась полная энергия системы. Параметры, по которым проверялась сходимость, подбирались таким образом, чтобы точность расчета полной энергии на атом составляла не менее, чем 1 мэВ. Эти потенцилы тестировались на сходимость путем запуска серии расчетов с разными наборами k-точек при нескольких значениях параметра Ecut в несколько этапов: 1. Сначала определялся оптимальный параметр Ecut при максимальной заданной плотности k-точек, затем, при выбранном Ecut , определялся оптимальный набор k-точек. 2. Полная оптимизация геометрии структуры методом Бройден-ФлетчерГольдфарб-Шано [63], при определенных оптимальных параметрах. 3. После оптимизации геометрии повторялся первый этап. Начальные значения параметра Ecut , для расчетов в программных пакете VASP составляли 1000 эВ, начальные значения плотности к-точек 0.08 Å. Заметим, что практика показала, как оптимальные значения параметра Ecut и наборы k-точек не менялись после релаксации параметров решетки структур чистых элементов. Наконец, полные релаксации геометрии 55

структур проводились до тех пор, пока силы, действующие на атомы не становились меньше 5 мэВ/Åпри расчетах в пакете VASP. В результате, параметры, по которым проверялась сходимость, подбирались таким образом, чтобы точность расчета полной энергии на атом составляла не хуже, чем 1 мэВ. Для расчета энергии образования межузельные атомов γ[111], γ[100] и γ[110] была построена свехъячейка алмаза состоящая из 216 атомов с данными видами дефектов. На рисунке 2.3, преставлен вид плоскости [110], красным цветом показаны смещеные атомы и межузельники. Рис. 2.3: Межузельные атомы в структуре алмаза, плоскость 110 Также была построена свехъячейки, для графена из 128 атомов, со смещениями 0.05 Å, 0.1 Å, 0.15 Å, в направлениях <0-10>, <001>, <010> и <210>, для графена. Используя программный пакет Phonopy были постороены смещения на 0.05 Å, 0.1 Å, 0.15 Å, для сверхъячейки алмаза, состоящей из 216 атомов. 56

2.2 Методика построение потенциалов межатомных взаимодействий Процесс построения потенциалов межатомных взаимодействий для углерода заключается в оптимизации значения потенциальных функций в узлах сплайнов относительно подгоночной базы данных, состоящей из экспериментальных величин и величин полученных вычислительным методом на основе теории функционала плотности. Потенциальные функции, представлятся в виде кубических сплайнов с равноудаленными узлами [], которые в свою очередь и являются оптимизируемыми параметрами функций. Для оптимизации функций Φ(R) на интервале [Rmin , RΦ ], использовали 30 узлов сплайна, для функций ρ(R) на интервале [Rmin , RΦ ] 25 узлов сплайна, при описании функций f3p было применено 15 узлов сплайна на интервале [Rmin , Rf ], в случаи функций g3pq 10 узлов сплайна на интервале [-1, +1] и 10 узлов сплайна для оптимизации функций F расположеной на интервале [0, 25]. Диапазон значений, для функции F(ρi ) включает все значения аргумента ρi , которые появляются только в процессе оптимизации. Для определения функций Φ, ρ, f3p был выбран внутренний радиус Rmin равный 1.1 Å. Внутренний радиус выбирается таким образом, чтобы он был меньше, чем ближайщее расстояние между атомами углерода в подгоночной базе данных, в нашем случае не меньше 1.43 Å-наименьшее растояние между атомами углерода в графене. За внешний радиус для функций парного потенциала и функций ρ было выбрано значение 4.4 Å. Внешний радиус для этих функций выбирался таким образом, чтобы он был больше, чем наибольшее расстояния между атомами углерода в подгоночной базе данных, в нашем случае больше 3.35 Åэто наибольшее растояние между атомами углерода в графите. Внешний радиус для базисных функций f3p , выбирался таким образом, что бы он был вдвое меньше внешнего радиуса парного потенциала и функций ρ и равен 2.2 Å. При этом целевую функцию Z мы минимизируем, используя метод минимизации, основанный на методах моделируемого отжига и Монте-Карло. Мы построили и протестировали четыре потенциала для углерода, с одной, двумя, тремя и четырьмя базисными функциями. В результате тестирования всех 57

построенных потенциалов, мы выбрем один, который далее будем обозначать, как POT C. Для POT C, чтобы точно учеть отталкивание на малых расстояних (r < RZBL ), мы используем отталкивающий потенциал ZieglerBiersack-Littmark (ZBL) [64], параметры для ηe и ηr в нашем случае равны 1,09863087 и 1,00648965, соответственно. На интервале RZBL < r < Rmin парные взаимодействия Φ(r) описываются полиномом пятой степени, коэффициенты которого обеспечивают непрерывность потенциала вместе с его первой и второй производной [54]. 2.3 Техника безопасности При работе с компьютером необходимо соблюдать следующие правила поведения: 1. Провода, розетки работающего компьютера и электрические вилки, руками не трогать. 2. При работе на компьютере необходимо, исключить попадание в него воды, поэтому не желательна работа в сырой одежде или с мокрыми руками. 3. При нарушении целостности корпуса комьютера необходимо прекратить работу. 4. Если появился запах гари или необычные звуки, необходимо немедленно выключить компьютер. 5. Нельзя оставлять компьютер без присмотра, поэтому при появлении в процессе работы, каких либо неотложных дел, необходимо выключить компьютер, если срок отсутствия превышает 20 мин . 6. Предметы на компьютер класть нежелательно, так как при этом уменьшается теплоотдача металлических элементов. Необходимо также соблюдать требования безопасности перед началом работы на компьютере. 1. Осмотреть ребочее место и в случае чего привести его порядок. 2. Необходимо отрегулировать освещение на рабочем месте и в убедится в том, что потока встречного света не будет препятсвовать работе. 3. Убедится в том, что электрооборудование правильно подключено к сети. 4. Поверхность экрана необходимо протереть салфеткой. 5. Установка стола и клавиатуры не должны приводить к некомфорту при работе. Последовательность включения компьютера: 1. Включить блок питания. 2. Включить периферийные устройства. 3. Включить системный блок. Требования безопасности во время работы. 1. Продолжительность работы перед 58

экраном не должна превышать 1 часа. 2. В течении всего рабочего времени стол содержать в порядке. 3. Открыть все вентиляционные устройства. 4. Выполнять санитарные нормы: соблюдать режим работы и отдыха. 5.Соблюдать правила эксплуатации и вычислительной техники в соответствии с инструкциями. 6. Соблюдать расстояние до экрана в пределах 70-80см. 7. Соблюдать установленный временем режим работы. Выполнять упражнения для рук, глаз и т.д. 8. Во время работы запрещается одновременно касаться экрана и клавиатуры. 9. Запрещается касаться задней панели системного блока при включённом питании. 10. Запрещается попадание воды на системный блок, рабочую поверхность и другие устройства. 11. Запрещается производить самостоятельное вскрытие и ремонт оборудования. 12 После работы на компьютере не рекомендуется смотреть телевизор 2-3часа. Требования безопасности в аварийных ситуациях 1. Во всех случаях обрывов проводов питания, неисправности заземления необходимо выключить компьютер. 2. В случае появления рези в глазах, резком ухудшении видимости, появлении боли в пальцах немедленно покинуть рабочее место сообщить руководителю работ и обратится к врачу. 3. При возгорании оборудования отключить питание и принять меры тушения. Требования безопасности после окончания работы 1. Произвести закрытие всех активных задач. 2. Выключить питание системного блока. 3. Выключить питание всех периферийных устройств 4. Отключить блок питания. 5. По окончанию работы осмотреть и привести в рабочее состояние компьютер. 59

Глава 3 Результаты исследования и их обсуждение 3.1 Тестирование PAW потенциалов График зависимости свободной энергии от объема ячейки, состоящей из восьми атомов, для структуры алмаза, представлен на рисунке 3.1. График имеет четкий минимум, энергия экспоненциально возрастает с уменьшением и увеличеним объема относительно этого зачения и хорошо согласует с уравнением состояния Vinet. Разность между экспериметальными результатами и данными полученными PAW потенциалами равновесного параметра решетки алмаза, составила 0.006 Å, -0.034 Å, 0.004 Å, 0.004 Å, 0.003 Å, -0.003 Å, для соотвествующих C, Сd, Ch, Cs, Chnr, CGW потенциалов. При этом различие по модулю всестороннего сжатия, отличется на -9.618 ГПа, 25.819 ГПа, -9.337 ГПа, 6.331 ГПа, -9.084 ГПа, 9.682 ГПа, соотвествнно, тем же C, Сd, Ch, Cs, Chnr, CGW потенциалам. Все потенциалы описывают параметр решетки алмаза в пределах ошибки эксперимента, наибольшее откланение имеет потенциал Cd. Также потенциал Cd дает наибольшую ошибку при описании модуля всесторонненого сжатия алмаза, наименьшей ошибкой обладает потенциал Cs. График зависимости свободной энергии от объема ячейки, состоящей из двух атомов, для структуры графена, представлена на рисунке 3.2. График 60

Рис. 3.1: Зависимость энергии от объема ячейки состоящей из восьми атомов для структуры Алмаза. Сплошная линия соответствует уравнению состояния Vinet. Таблица 3.1: Расчет свойств аламаза и графена, потенциалами ряда C, Сd, Ch, Cs, Chnr, CGW Решетка Cвойство С Cd Ch Cs Chn r CGW Эксп. Алмаз 3.53 3.57 3.57 3.57 Графен a0 (Å) 3.57 B(ГПа) 435.38 470.82 435.66 438.67 435.92 435.32 442 a0 (Å) 2.47 2.45 2.76 2.47 2.47 3.57 2.47 3.57 2.46 имеет четкий минимум, энергия экспоненциально возрастает с уменьшением и увеличеним объема относительно этого значения и хорошо согласуется с уравнением состояния Vinet. Данные рассчетов равновесного параметра решетки и модуля всестороннего сжатия для алмаза и равновесный параметр решетки для графена, полученого из PAW потенциала, представленны в таблице 3.1. Разницы энергий между структурами алмаза-графит и алмаз-графен представлены в таблице 3.2 и имеют следующие значения для ряда потенциалов C, Сd, Ch, Cs, Chnr, CGW: 0.123 эВ, -0.010 эВ, 0.122 эВ, 0.109 эВ, 0.121 эВ, 0.119 эВ и 0.137 эВ, -0.034 эВ, 0.135 эВ, 0.124 эВ, 0.134 эВ, 0.133 61

Рис. 3.2: Кривая зависимости общей энергии от объема, для графена Таблица 3.2: Расчет свойств аламаза и графена, потенциалами ряда C, Сd, Ch, Cs, Chn, CGW δE(эВ) С Cd Ch Cs Chn r CGW Графит 0.123 -0.010 0.122 0.109 0.121 0.119 Графен 0.137 -0.034 0.135 0.124 0.134 0.133 эВ. Алмаз является метастабильной фазой и должен обладать большей энергией сублимации, по отношению к стабильным фазам графиту и графену, чему и соотвествуют расчеты разницы энергий сублимаций этих структур потенциалами C, Ch, Cs, Chnr и CGW. Потенциал Cd показывает неверный результ, представляя алмаз, как более стабильную структуру, чем графит и графен. Возможно это связано с наличием большего количества плоских волн в нем, чем в других потенциалах. Одной из важнейших характеристик при выборе PAW потенциалов, является скорость расчетов. На рисунке 3.3 предстален график зависимости скорости расчета ячейки алмаза состоящей из восьми атомов, от используемого PAW потенциала Наибольшей скоростью расчета обладает PAW потенциал с на62

Рис. 3.3: Скорости расчета ячейки алмаза состоящей из восьми атомов, взависимости от используемого PAW потенциала званием Cd, далее идут в порядке убывания скорости расчета потенциалы Ch, C, CGW, Cs, Chnr. Исходя из представленных данных, мы считаем, что оптипальным, для дальнейшего исследования является PAW потенциал С. С целью уменьшения времени расчета, была проведена сходимость по k точкам и сходимость по ENCAT. В таблице 3.3 представлена разность энергий алмаз-графен взависимости от плотности к-точек. Из представленых данных мы видим, что при 0.2 Å−1 к-точек, энергия начинает падать, в связи с чем была выбрана эта плотность, для дальнейших расчетов. В таблице 3.4 представлена разность энергий алмаз-графен взависимости от энергии обрезания, с заданой плотностью к-точек, выбранной ранее(0.2 Å−1 ). Из представленых данных мы видим, что при 300 эВ, разность энергий алмаз-графен начинает падать. Для больше точности расчетов была выбрана энергия обрезания в 400 эВ и использовалась, для дальнейших расчетов. В результате проверки сходимости по k- точкам и энергии обрезания, k0.2 был выбран PAW потенциал CEN 400 , c энергией обрезания 400 эВ и плот- ностью k-точек 0,2 Å−1 График зависимости свободной энергии от объема ячейки, состоящей из шестидесяти атомов молекулы фуллерена, представлен на рисунке 3.4. График имеет четкий минимум, энергия экспоненциально возрастает с умень63

Таблица 3.3: Расчет свойств аламаза и графена, потенциалами ряда C, Сd, Ch, Cs, Chn, CGW KSPACING, (Å−1 ) E0 [алмаз]-E0 [графен], (эВ) 0.08 0.132 0.12 0.134 0.10 0.133 0.14 0.134 0.16 0.134 0.18 0.134 0.20 0.132 0.22 0.134 0.24 0.134 Таблица 3.4: Сходимость разности энергии алмаз-графен взависимости от энергии обрезания ENCAT, (эВ) E0 [алмаз]-E0 [графен], (эВ) 1000 0.132 900 0.131 800 0.131 700 0.133 600 0.130 500 0.129 400 0.133 300 0.103 64

Рис. 3.4: Расчет энергии сублимации, для фуллерена шением и увеличеним объема относительно этого значения и хорошо согласуется с уравнением состояния Vinet. Энергия сублимации фуллерена отличается от алмаза на 0.24 эВ, в пользу стабильности алмаза. Диаметр молекулы фуллерена составил 3.75 Å, при этом отличаясь от экспериментальных значений на 0.05 Å. 3.2 Подгоночная база данных, для проверки на сходимость характеристик структур углерода В таблице 3.5 представлены параметры решётки и атомные разности энергий относительно равновесной решётки алмаза для структур из подгоночной базы данных, рассчитаные в рамкх ТФП. 65

Таблица 3.5: Параметры решётки и атомные разности энергий относительно равновесной решётки алмаза для структур из подгоночной базы данных Решетка Метод a0 c0 V0 ∆E ОЦК GGA 2.17 2.17 5.14 4.207 ГЦК GGA 2.79 2.79 6.55 4.422 ПК GGA 1.77 1.77 5.59 2.622 Графен GGA 2.47 20 0.109 0.133 Эксп. 2.47 - - - GGA 2.46 6.70 8.78 0.019 Эксп. 2.46 6.70 8.78 0.004 Наш GGA 3.57 3.57 5.70 Эксп. 3.57 5.70 Графит Алмаз 3.57 3.3 Расчет сходимости в описании свойств алмаза и графита На рисунке 3.5 представлены потенциальные функции, описывающие межатомные взаимодействия, для углерода с одной базисной функцией. В этом случае, межатомные взаимодействия, описываются пятью потенциальными функциями: Φ(R), f31 (R), g311 (cos(θ)), ρ(R) , F (ρ). 66

Рис. 3.5: Потенциальные функции Φ(r) (в эВ), ρ(r) (безразмерная), F (ρ̄) (в эВ), f p (r) (безразмерные) и g pq (x) (в эВ) взаимодействий между атомами углерода, задаваемые кубическими сплайнами. Значения функций в узлах сплайнов (параметры потенциальных функций) обозначены символами. На рисунке 3.6 представлены потенциальные функции, описывающие межатомные взаимодействия для углерода с двумя базисными функциями. В этом случае, межатомные взаимодействия, описываются восьмью потенциальными функциями: Φ(R), f31 (R), f32 (R), g311 (cos(θ)), g312 (cos(θ)), g322 (cos(θ)), g333 (cos(θ)), ρ(R) , F (ρ). На рисунке 3.7 представлены потенциальные функции, описывающие межатомные взаимодействия для углерода с тремя базисными функциями. В этом случае, межатомные взаимодействия, описываются двеннадцатью потенциальными функциями: Φ(R), f31 (R), f32 (R), f33 (R), g311 (cos(θ)), g312 (cos(θ)), g322 (cos(θ)), g333 (cos(θ)), g334 (cos(θ)), g344 (cos(θ)), ρ(R) , F (ρ). На рисунке 3.8 представлены потенциальные функции, описывающие межатомные взаимодействия для углерода с четырьмя базисными функциями. В этом случае, межатомные взаимодействия, описываются семнадцатью потенциальными функциями: Φ(R), f31 (R), f32 (R), f33 (R), f34 (R), g 1 13 (cos(θ)), g312 (cos(θ)), g322 (cos(θ)), g333 (cos(θ)), g334 (cos(θ)), g344 (cos(θ)), g345 (cos(θ)), 67

Рис. 3.6: Потенциальные функции Φ(r) (в эВ), ρ(r) (безразмерная), F (ρ̄) (в эВ), f p (r) (безразмерные) и g pq (x) (в эВ) взаимодействий между атомами углерода, задаваемые кубическими сплайнами. Значения функций в узлах сплайнов (параметры потенциальных функций) обозначены символами. g355 (cos(θ)), g356 (cos(θ)), g366 (cos(θ)), ρ(R) , F (ρ). На представленных выше рисунках 3.5, 3.6, 3.7, 3.8, для функций описывающих парное взаимодействие Φ(r) наблюдается увеличение значения функции с уменьшением расстояние между атомами, что хорошо согласуется с экспериментальными данными, при описании данного вида взаимодействия. Функции F(ρ), имеют по одному минимуму и обладают линейной зависимостью, в случаях когда аргумент принимает большие значения и стремится к постоянным значениям. Полученный вид фукнции подобен функциям погружения в теории эффективной среды [], которая была использована для обоснования МПА []. Так же следует заметить, что функция ρ(r), может принимать как положительные, так и отрицательные значения, в отличие от функции электронной плотности в МПА, которая принимает только положительные значения. На рисунке 3.9 представлена зависимость модулей упругости алмаза, от числа базисных функций и их сравнение с экспериментальными данными. 68

Рис. 3.7: Потенциальные функции Φ(r) (в эВ), ρ(r) (безразмерная), F (ρ̄) (в эВ), f p (r) (безразмерные) и g pq (x) (в эВ) взаимодействий между атомами углерода, задаваемые кубическими сплайнами. Значения функций в узлах сплайнов (параметры потенциальных функций) обозначены символами. На рисунке 3.10 представлена зависимость модулей упругости алмаза, от числа базисных функций и их сравнение с экспериментальными данными. В результате проверки на сходимость потенциала с одной базисными функциями в трехчатичных взаимодействиях, при описании констант упругости в структурах алмаза и графита, константы упругости, для графита C11 , C12 , C13 , C33 и C66 , отличаются на 174.86 ГПа, 195.41 ГПа, 15.60 ГПа, 21.76 ГПа, 0.17 ГПа и 58.13 ГПа, соотвественно, константы упругости, для алмаза C11 , C12 и C44 отличаются на 38.35 ГПа, 4.27 ГПа, 0 ГПа, соотвественно, от их экспериментальных значений. В результате проверки на сходимость потенциала с двумя базисными функциями в трехчатичных взаимодействиях, при описании констант упругости в структурах алмаза и графита, константы упругости, для графита C11 , C12 , C13 , C33 , C44 и C66 , отличаются на 57.70 ГПа, 6.57 ГПа, 3.75 ГПа, 3.81 ГПа, 0.11 ГПа и 33.13 ГПа, соотвественно, константы упругости, 69

Рис. 3.8: Потенциальные функции Φ(r) (в эВ), ρ(r) (безразмерная), F (ρ̄) (в эВ), f p (r) (безразмерные) и g pq (x) (в эВ) взаимодействий между атомами углерода, задаваемые кубическими сплайнами. Значения функций в узлах сплайнов (параметры потенциальных функций) обозначены символами. для алмаза C11 , C12 и C44 отличаются на 50.26 ГПа, 5.01 ГПа, 16.03 ГПа, соотвественно, от их экспериментальных значений. В результате проверки на сходимость потенциала с тремя базисными функциями в трехчатичных взаимодействиях, при описании констант упругости в структурах алмаза и графита, константы упругости, для графита C11 , C12 , C13 , C33 , C44 и C66 , отличаются на 33.00 ГПа, 3.75 ГПа, 0.09 ГПа, 2.03 ГПа, 0.08 ГПа и 17.94 ГПа, соотвественно, константы упругости, для алмаза C11 , C12 и C44 отличаются на 47.79 ГПа, 3.75 ГПа, 77.16 ГПа, соотвественно, от их экспериментальных значений. В результате проверки на сходимость потенциала с четырьмя базисными функциями в трехчатичных взаимодействиях, при описании констант упругости в структурах алмаза и графита, константы упругости, для графита C11 , C12 , C13 , C33 , C44 и C66 , отличаются на 23.17 ГПа, 0.48 ГПа, 3.22 ГПа, 0.53 ГПа, 0.26 ГПа и 11.49 ГПа, соотвественно, константы упругости, для алмаза C11 , C12 и C44 отличаются на 10.81 ГПа, 6.38 ГПа, 13.75 ГПа, 70

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

Рис. 3.10: Расчет констант упругости с использованием межатомного потенциала, для графита На рисунке 3.11 представлен график зависимости модуля разности |CijP OT CijExp | констант упругости графита, от числа базисных функций. На рисунке 3.12 представлен график зависимости модуля разности |CijP OT CijExp | констант упругости алмаза, от числа базисных функций. 72

Рис. 3.11: Зависимости модуля разности |CijP OT -CijExp | констант упругости графита, от числа базисных функций Рис. 3.12: Зависимости модуля разности |CijP OT -CijExp | констант упругости алмаза, от числа базисных функций 73

Рис. 3.13: Скорость расчета параметра сверхъячейки алмаза состоящей из 2744 атомов взависимости от числа базисных функций. Время расчета параметра сверхъячейки алмаза в 2744 атомов, для потенциалов с одной, двумя, тремя и четырьмя базисными фукнциями составило 12.02 с, 13.02 с, 16.42 с, 23.07 с, соотвественно. На рисунке 3.13 представлен график зависимости скорости расчета параметра сверхъячейки алмаза состоящей из 2744 атомов, от числа базисных функций. Набюлюдается снижение вычислительной скорости с увеличеним числа базисных функций, причем при двух базисных функциях скорость расчета уменьшается в 1.08 раза, с тремя в 1.37 раза и с четырьмя в 1.92 раза, относительно потенциала с одной базисной функцией. Из данного исследование на сходимость, по модулям упругости алмаза и графита, учитывая, также время расчета, целесообразней использовать три базисных функций. Поэтому дальнейшие вычисления были проведены с потенциалом имеющим три базисные функции. 74

3.4 Расширенная подгоночная база данных, для потенциала В таблице 3.6 представлены параметры решётки и атомные разности энергий относительно равновесной решётки алмаза для структур из подгоночной базы данных. Таблица 3.6: Параметры решётки и атомные разности энергий относительно равновесной решётки алмаза для структур из подгоночной базы данных Решетка Метод a0 c0 V0 ∆E ОЦ-8 GGA 4.43 4.43 5.44 0.851 α-кб GGA 2.57 20 514 0.8727 β-кб GGA 2.57 20 514 0.8716 Лонсдейлит GGA 2.52 4.18 - 0.0196 Эксп. 2.41 4.17 5.30 - В таблице 3.7 представлены экспериментальные данные по равновесной решётке графита и константам упругости, использованные для оптимизации параметров потенциалов. 75

Таблица 3.7: Экспериментальные данные по равновесной решётке графита и константам упругости, использованные для оптимизации параметров потенциалов. Ec - энергия сублимации, a0 - параметр решётки, Cij - модули упругости Свойство Значение Ec , эВ 7.374 a0 , Å 2.46 c/a 2.72 C11 , ГПа 1109 C12 , ГПа 139 C13 , ГПа 0 C33 , ГПа 38 C44 , ГПа 5 C66 , ГПа 485 На рисунке 3.14 Расчет равновесных параметров решетки , для Лонсдейлита в рамках теории функционала плотности. Отметим, что диапазон значений энергий, для парметров решетки лежит в центре рисунка. 76

Рис. 3.14: Расчет равновесных параметров решетки , для Лонсдейлита 3.5 Тестирование межатомного потенциала На рисунке 3.15 представлены потенциальные функции, описывающие межатомные взаимодействия для углерода с тремя базисными функциями. В этом случае, межатомные взаимодействия, описываются двеннадцатью потенциальными функциями: Φ(R), f31 (R), f32 (R), f33 (R), g311 (cos(θ)), g312 (cos(θ)), g322 (cos(θ)), g333 (cos(θ)), g334 (cos(θ)), g344 (cos(θ)), ρ(R) , F (ρ). На представленном выше рисунке 3.15, для функции описывающих парное взаимодействие Φ(r) наблюдается увеличение значения функции с уменьшением расстояние между атомами, что хорошо согласуется с экспериментальными данными, при описании данного вида взаимодействия. Функции F(ρ), имеют по одному минимуму и обладают линейной зависимостью, в случаях когда аргумент принимает большие значения и стремится к постоянным значениям. Полученный вид фукнции подобен функциям погружения в теории эффективной среды, которая была использована для обоснования МПА. Так же следует заметить, что функция ρ(r), может принимать как положительные, так и отрицательные значения, в отличие от функции электронной плотности в МПА, которая принимает только положительные значения. В таблице 3.8 представленны данные об константах упругости для алма- 77

Рис. 3.15: Потенциальные функции Φ(r) (в эВ), ρ(r) (безразмерная), F (ρ̄) (в эВ), f p (r) (безразмерные) и g pq (x) (в эВ) взаимодействий между атомами углерода, задаваемые кубическими сплайнами. Значения функций в узлах сплайнов (параметры потенциальных функций) обозначены символами. 78

за и графита, рассчитаные потенциалами MEAM, Терсоффа, REBO, BOP, ALREBO и POT C в сравнении с экспериментальными данными. Потенциалы MEAM и Терсоффа дают ошибку при описании константы упругости С12 в грифите в виде отрицательного значения. Но следует отметить, что MEAM лучше других потенциалов описывает константы упругости, для алмаза в сравнении с экспериментом, отличаясь на 1 ГПа, 2 ГПа и 44 ГПа, для С11, С12 и С44, соответственно. Таблица 3.8: Упругие характеристики, для алмаза и графена, рассчитаные потенциалами MEAM, Терсоффа, REBO, BOP и POT C, в сравнении с экспериментальными данными. Crystal Diamond Graphite Method C11 C12 C13 C33 C44 C66 B Exp. 1080 125 579 442 MEAM[10] 1079 127 623 444 BOP[40] 1073 224 609 507 REBO [65] 1070 120 720 - ALREBO [37] 1120 130 720 - Терсофф [66] 1090 120 640 443 POT C 1065 128 578 440 Exp. 1109 139 0 MEAM+LJ[10] 1099 -450 0.03 BOP[40] 38 5 485 38 0.03 - 961 497 - - - 240 REBO[65] 1346 497 - - - - ALREBO[37] 1150 150 10 40 - - Терсофф [5] 1210 -190 - - - - 139 0.35 39 3.8 489 POT C 1116 36 37 В таблице 3.9 представлены упругие свойства графита и алмаза, рассчитанные с использованием построенных потенциалов, в сопоставлении с экспериментальными данными. Энергия образования дефекта γ[111] отличается от экспериментальных значений на 3.33 эВ и на 1.08 эВ от ТФП расчета, при этом обладает наименьшей энергий относительно других де79

фектов, представленных в таблице 3.9, что хорошо согласуется с экспериментальными и ТФП данными. Энергия образования дефекта γ[110], отличается от экспериментальных значений на 4.76 эВ и на 1.51 эВ от ТФП расчета. Энергия образования дефекта γ[100] отличается от экспериментальных значений на 2.16 эВ и на 1.03 эВ от ТФП расчета, при этом обладает наибольшей энергий относительно других дефектов, представленных в таблице 3.9, что хорошо согласуется с экспериментальными и ТФП данными. Таблица 3.9: Энергия образования дефектов в алмазе Метод Тип дефекта ПОТ. γ[111] 8.63 γ[110] 11.36 γ[100] 11.36 γ[111] 5.3 γ[110] 6.5 γ[100] 9.2 γ[111] 9.71 γ[110] 9.77 γ[100] 10.33 Эксп.[67] ТФП. Энергия образования, эВ График зависимости свободной энергии от объема ячейки, состоящей из шестидесяти атомов молекулы фуллерена, представлен на рисунке 3.16. График имеет четкий минимум, энергия экспоненциально возрастает с уменьшением и увеличеним объема относительно этого значения и хорошо согласуется с уравнением состояния Vinet. Энергия сублимации фуллерена отличается от алмаза на 0.3 эВ, в пользу стабильности алмаза. Отличие в разницах энергии фуллерена относительно результатов ТФП, составляет 0.06 эВ и находится из следующего соотношения |∆E DF T − ∆E P OT |. Диаметр молекулы фуллерена составил 3.75 Å, при этом отличаясь от экспериментальных значений на 0.2 Åи на 0.18 Åот ТФП расчетов. 80

Рис. 3.16: Фуллерен С60 81

Заключение 1. В рамках теории функционала плотности, были рассчитаны энергии сублимации и параметров решеток структур алмаза, графита, графена, ОЦК и ГЦК, для подгоночной базы данных межатомных потенциалов. Полученные данные имеет хорошее согласие с исследованиями других авторов, что говорит об их надежности. 2. В рамках нового подхода были построены потенциалы, для углерода, с одной, двумя, тремя и четырьмя базисными функциями в трехчастичных взаимодействиях. Для оптимизации функций Φ(R) на интервале [Rmin , RΦ ], использовали 30 узлов сплайна, для функций ρ(R) на интервале [Rmin , RΦ ] 25 узлов сплайна, при описании функций f3p было применено 15 узлов сплайна на интервале [Rmin , Rf ], в случаи функций g3pq 10 узлов сплайна на интервале [-1, +1] и 10 узлов сплайна для оптимизации функций F расположеной на интервале [0, 25]. 3. В результате сравнительного анализа было выяснено, что необходимо использовать три базисных функции, при этом расчет констант упругости, для графита находится в пределах ошибки и не превышает 31 %, для алмаза не превышает 0.5%. Ранее построенный потенциал, основаный на методе погруженного атома MEAM, для углерода, дает ошибки при расчете констант упругости графита в 263%, а для алмаза в 7.5%. 4. При тестировании межатомного потенциала с тремя базисными функциями, было выяснено, что полученный потенциал, находится в хорошем согласии с данными полученными в рамках теории функци82

онала плотности, при описании точечных дефектов в алмазе и дает меньшую ошибку 16%, в отличии от потенциала типа MEAM, ошибка которого составляет 30%. Также были расчитаны энергии сублимации и диаметр фуллерена, значения которых составили 7.13 эВ и 3.75 Å, соотвественно. 83

Литература [1] Dynamics of radiation damage / J. B. Gibson, A. N. Goland, M. Milgram, G. H. Vineyard // Physical Review. — 1960. — Vol. 120, no. 4. — Pp. 1229– 1253. [2] Daw Murray S., Baskes M. I. Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals // Physical Review Letters. — 1983. — Vol. 50, no. 17. — Pp. 1285–1288. [3] Baskes M. I. Application of the Embedded-Atom Method to Covalent Materials: A Semiempirical Potential for Silicon // Physical Review Letters. — 1987. — Vol. 59, no. 23. — Pp. 2666–2669. [4] Stillinger Frank H., Weber Thomas A. Computer simulation of local order in condensed phases of silicon // Physical Review B. — 1985. — Vol. 31, no. 8. — Pp. 5262–5271. [5] Tersoff J. Empirical interatomic potential for carbon, with applications to amorphous carbon // Physical Review Letters. — 1988. — Vol. 61, no. 25. — Pp. 2879–2882. [6] A Second-Generation Reactive Empirical Bond Order (Rebo) Potential Energy Expression for Hydrocarbons / D W Brenner, O A Shenderova, J A Harrison et al. // J. Phys.: Condens. Matter. — 2002. — Vol. 14. — P. 783. [7] Simulating radiation damage cascades in graphite / H. J. Christie, M. Robinson, D. L. Roach et al. // Carbon. — 2015. — Vol. 81, no. 1. — Pp. 105–114. 84

[8] Molecular dynamics simulation of radiation damage cascades in diamond / J. T. Buchan, M. Robinson, H. J. Christie et al. // Journal of Applied Physics. — 2015. — Vol. 117, no. 24. [9] Anomalous Dynamical Behavior of Freestanding Graphene Membranes / M. L. Ackerman, P. Kumar, M. Neek-Amal et al. // Physical Review Letters. — 2016. — Vol. 117, no. 12. — Pp. 1–5. [10] Lee Byeong Joo, Lee Jin Wook. A modified embedded atom method interatomic potential for carbon // Calphad: Computer Coupling of Phase Diagrams and Thermochemistry. — 2005. — Vol. 29, no. 1. — Pp. 7–16. [11] Lipnitskii A. G., Saveliev V. N. Development of n-body expansion interatomic potentials and its application for V // Computational Materials Science. — 2016. — Vol. 121. — Pp. 67–78. [12] Kohn W., Sham L. J. Self-Consistent Equations Including Exchange and Correlation Effects // Physical Review. — 1965. — 11. — Vol. 140, no. 4A. — Pp. A1133–A1138. — URL: http://link.aps.org/doi/10.1103/ PhysRev.140.A1133. [13] Martin R.M. Electronic Structure: Basic Theory and Practical Methods. — Cambridge: Cambridge University Press, 2004. [14] Perdew John P., Wang Yue. Accurate and simple analytic representation of the electron-gas correlation energy // Physical Review B. — 1992. — 6. — Vol. 45, no. 23. — Pp. 13244–13249. — URL: http://link.aps.org/ doi/10.1103/PhysRevB.45.13244. [15] Perdew John P., Burke Kieron, Ernzerhof Matthias. Generalized Gradient Approximation Made Simple // Physical Review Letters. — 1996. — 10. — Vol. 77, no. 18. — Pp. 3865–3868. — URL: http://link.aps.org/ doi/10.1103/PhysRevLett.77.3865. [16] Blöchl P. E. Projector augmented-wave method // Physical Review B. — 1994. — 12. — Vol. 50, no. 24. — Pp. 17953–17979. — URL: http: //link.aps.org/doi/10.1103/PhysRevB.50.17953. 85

[17] Alder B. J., Wainwright T. E. Studies in molecular dynamics. I. General method // The Journal of Chemical Physics. — 1959. — Vol. 31, no. 2. — Pp. 459–466. [18] Verlet Loup. Computer "Experiments"on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules // Physical Review. — 1967. — 7. — Vol. 159, no. 1. — Pp. 98–103. — URL: http://link.aps.org/doi/10.1103/PhysRev.159.98. [19] Jones J. E. On the Determination of Molecular Fields. II. From the Equation of State of a Gas // Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. — 1924. — Vol. 106, no. 738. — Pp. 463–477. — URL: http://rspa.royalsocietypublishing. org/cgi/doi/10.1098/rspa.1924.0082. [20] Buckingham R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon // Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. — 1938. — Vol. 168, no. 933. — Pp. 264– 283. — URL: http://rspa.royalsocietypublishing.org/cgi/doi/ 10.1098/rspa.1938.0173. [21] Jones J. E. On the Determination of Molecular Fields. II. From the Equation of State of a Gas // Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. — — Vol. 106, no. 738. — Pp. 463–477. — 1924. — 10. URL: http://rspa. royalsocietypublishing.org/cgi/doi/10.1098/rspa.1924.0082. [22] Tersoff J. New empirical model for the structural properties of silicon // Physical Review Letters. — 1986. — 2. — Vol. 56, no. 6. — Pp. 632–635. — URL: http://link.aps.org/doi/10.1103/PhysRevLett.56.632. [23] Elastic constants of compression-annealed pyrolytic graphite / O. L. Blakslee, D. G. Proctor, E. J. Seldin et al. // Journal of Applied Physics. — 1970. — Vol. 41, no. 8. — Pp. 3373–3382. [24] Marks N. A. Generalizing the environment-dependent interaction potential 86

for carbon // Physical Review B - Condensed Matter and Materials Physics. — 2001. — Vol. 63, no. 3. — Pp. 1–7. [25] Adelman S. A. Generalized Langevin equation approach for atom/solidsurface scattering: Collinear atom/harmonic chain model // The Journal of Chemical Physics. — 1974. — Vol. 61, no. 10. — P. 4242. — URL: http://scitation.aip.org/content/aip/journal/jcp/61/10/ 10.1063/1.1681723. [26] Andersen Hans C. Molecular dynamics simulations at constant pressure and/or temperature // The Journal of Chemical Physics. — 1980. — Vol. 72, no. 4. — P. 2384. — URL: http://scitation.aip.org/content/ aip/journal/jcp/72/4/10.1063/1.439486. [27] Molecular dynamics with coupling to an external bath / H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren et al. // The Journal of Chemical Physics. — 1984. — Vol. 81, no. 8. — P. 3684. — URL: http://scitation.aip.org/content/aip/journal/jcp/81/8/ 10.1063/1.448118. [28] Nose Shuichi. A molecular dynamics method for simulations in the canonical ensemble // Molecular Physics. — 1984. — 8. — Vol. 52, no. 2. — Pp. 255–268. — URL: http://www.tandfonline.com/doi/abs/10. 1080/00268978400101201. [29] Tersoff J. Modeling solid-state chemistry: Interatomic potentials for multicomponent systems // Physical Review B. — 1989. — Vol. 39, no. 8. — Pp. 5566–5568. [30] Tersoff J. Chemical order in amorphous silicon carbide // Physical Review B. — 1994. — Vol. 49, no. 23. — Pp. 16349–16352. [31] Tang Meijie, Yip Sidney. Atomistic simulation of thermomechanical properties of β-SiC // Physical Review B. — 1995. — Vol. 52, no. 21. — Pp. 15150–15159. 87

[32] Nordlund K., Keinonen J., Mattila T. Formation of ion irradiation induced small-scale defects on graphite surfaces // Physical Review Letters. — 1996. — Vol. 77, no. 4. — Pp. 699–702. [33] Titantah J. T., Lamoen D. Sp3/sp2 characterization of carbon materials from first-principles calculations: X-ray photoelectron versus high energy electron energy-loss spectroscopy techniques // Carbon. — 2005. — Vol. 43, no. 6. — Pp. 1311–1316. [34] Brenner Dw. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films // Physical Review B. — 1990. — Vol. 42, no. 15. — Pp. 9458–9471. — URL: http://prb.aps. org/abstract/PRB/v42/i15/p9458_1. [35] Stuart Steven J., Tutein Alan B., Harrison Judith A. A reactive potential for hydrocarbons with intermolecular interactions // Journal of Chemical Physics. — 2000. — Vol. 112, no. 14. — Pp. 6472–6486. [36] Screened empirical bond-order potentials for Si-C / Lars Pastewka, Andreas Klemenz, Peter Gumbsch, Michael Moseler // Physical Review B - Condensed Matter and Materials Physics. — 2013. — Vol. 87, no. 20. — Pp. 1–12. [37] O’Connor Thomas C., Andzelm Jan, Robbins Mark O. AIREBO-M: A reactive model for hydrocarbons at extreme pressures // Journal of Chemical Physics. — 2015. — Vol. 142, no. 2. [38] Pettifor D. G., Oleinik I. I. Analytic bond-order potentials beyond tersoffbrenner. i. theory // Physical Review B - Condensed Matter and Materials Physics. — 1999. — Vol. 59, no. 13. — Pp. 8487–8499. [39] Atomistic modeling of hydrocarbon systems using analytic bond-order potentials / M. Mrovec, M. Moseler, C. Elsässer, P. Gumbsch // Progress in Materials Science. — 2007. — Vol. 52, no. 2-3. — Pp. 230–254. [40] Zhou X W, Ward D K, Foster M E. An Analytical Bond-Order Potential for Carbon. — 2015. — Pp. 1719–1735. 88

[41] ReaxFF: A reactive force field for hydrocarbons / a. C T Van Duin, Siddharth Dasgupta, Francois Lorant, William a. Goddard // J. Phys. Chem. A. — 2001. — Vol. 105. — Pp. 9396–9409. [42] Reaction Charging Chemical. Harding Battery Handbook. — Norton Shores. — Pp. 10–26. — URL: http://hardingenergy.com/handbook/. [43] Misra A, Demkowicz MJ. The radiation damage tolerance of ultra-high strength nanolayered composites // JOM Journal of the . . . . — 2007. — no. September. — Pp. 3–6. — URL: http://www.springerlink.com/ index/n566256ttw573m40.pdf. [44] Murrell J.N., Carter S., Farantos S.C. Molecular Potential Energy Function. — Chichester: John Wiley & Sons, Inc., 1984. — P. 206. [45] Many-body forces and electron correlation in small metal clusters / Ilya G. Kaplan, Jorge Hernández-Cobos, Octavio Novaro // Physical Review A. — Iván Ortega-Blake, 1996. — 4. — Vol. 53, no. 4. — Pp. 2493–2500. — URL: http://link.aps.org/doi/10.1103/ PhysRevA.53.2493. [46] Daw Murray S., Baskes M. I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals // Physical Review B. — 1984. — 6. — Vol. 29, no. 12. — Pp. 6443–6453. — URL: http://link.aps.org/doi/10.1103/PhysRevB.29.6443. [47] Baskes M. I. Many-Body Effects in fcc Metals: A Lennard-Jones Embedded-Atom Potential // Physical Review Letters. — 1999. — 9. — Vol. 83, no. 13. — Pp. 2592–2595. — URL: http://link.aps.org/doi/ 10.1103/PhysRevLett.83.2592. [48] Chantasiriwan Somchart, Milstein Frederick. Embedded-atom models of 12 cubic metals incorporating second- and third-order elastic-moduli data // Physical Review B. — 1998. — 9. — Vol. 58, no. 10. — Pp. 5996–6005. — URL: http://link.aps.org/doi/10.1103/PhysRevB.58.5996. 89

[49] Daw Murray S., Baskes M. I. Semiempirical, Quantum Mechanical Calculation of Hydrogen Embrittlement in Metals // Physical Review Letters. — 1983. — 4. — Vol. 50, no. 17. — Pp. 1285–1288. — URL: http://link.aps.org/doi/10.1103/PhysRevLett.50.1285. [50] Daw Murray S., Foiles Stephen M., Baskes Michael I. The embedded-atom method: a review of theory and applications // Materials Science Reports. — 1993. — 3. — Vol. 9, no. 7-8. — Pp. 251–310. — URL: http:// linkinghub.elsevier.com/retrieve/pii/092023079390001U. [51] Hu Wangyu, Shu Xiaolin, Zhang Bangwei. Point-defect properties in body-centered cubic transition metals with analytic EAM interatomic potentials // Computational Materials Science. — 2002. — 4. — Vol. 23, no. 1-4. — Pp. 175–189. — URL: http://linkinghub.elsevier.com/ retrieve/pii/S0927025601002385. [52] Stillinger Frank H., Weber Thomas A. Computer simulation of local order in condensed phases of silicon // Physical Review B. — 1985. — 4. — Vol. 31, no. 8. — Pp. 5262–5271. — URL: http://prb.aps.org/abstract/PRB/v31/i8/p5262{_}1http:// link.aps.org/doi/10.1103/PhysRevB.31.5262. [53] Mistriotis AD, Flytzanis N, Farantos SC. Potential model for silicon clusters // Physical Review B. — 1989. — 1. — Vol. 39, no. 2. — Pp. 1212–1218. — URL: http://prb.aps.org/abstract/PRB/v39/i2/ p1212{_}1http://link.aps.org/doi/10.1103/PhysRevB.39.1212. [54] Lipnitskii A.G., Saveliev V.N. Development of n-body expansion interatomic potentials and its application for V // Computational Materials Science. — 2016. — 8. — Vol. 121. — Pp. 67–78. — URL: http: //linkinghub.elsevier.com/retrieve/pii/S0927025616301549. [55] Baskes M I, Srinivasan S G. The embedded atom method ansatz: validation and violation // Modelling and Simulation in Materials Science and Engineering. — 2014. — 3. — 90 Vol. 22, no. 2. — P. 025025.

— URL: http://stacks.iop.org/0965-0393/22/i=2/a=025025?key= crossref.466d3328a9a1029ec0f1beab0f07dd37. [56] Kirkpatrick S., Gelatt C. D., Vecchi M. P. Optimization by Simulated Annealing // Science. — 1983. — 5. — Vol. 220, no. 4598. — Pp. 671–680. — URL: http://www.sciencemag.org/cgi/doi/10.1126/ science.220.4598.671. [57] Ab initio study of dihydrogen binding in metal-decorated polyacetylene for hydrogen storage / Hoonkyung Lee, Woon Choi, Manh Nguyen et al. // Physical Review B. — 2007. — 11. — Vol. 76, no. 19. — Pp. 1–7. — URL: http://link.aps.org/doi/10.1103/PhysRevB.76.195110. [58] Bottka N. 15 APRIL 1975 Theory. — 1975. — Vol. 11, no. 8. [59] Kittel Charles. Introduction to solid state physics / Charles Kittel. — Singapore : Wiley, 1986. — P. 646. — URL: http://trove.nla.gov. au/work/9250395. [60] Thermal expansion of a high purity synthetic diamond single crystal at low temperatures / Toshimaro Sato, Kazutoshi Ohashi, Tomoko Sudoh et al. // Physical Review B. — 2002. — Vol. 65, no. 9. — P. 092102. — URL: https://link.aps.org/doi/10.1103/PhysRevB.65.092102. [61] Kresse G. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set // Physical Review B. — 1996. — 10. — Vol. 54, no. 16. — Pp. 11169–11186. — URL: http://link.aps.org/ doi/10.1103/PhysRevB.54.11169. [62] Methfessel M., Paxton A. T. High-precision sampling for Brillouin-zone integration in metals // Physical Review B. — 1989. — 8. — Vol. 40, no. 6. — Pp. 3616–3621. — URL: http://link.aps.org/doi/10.1103/ PhysRevB.40.3616. [63] Head John D., Zerner Michael C. A Broyden—Fletcher—Goldfarb—Shanno optimization procedure for molecular geometries // Chemical Physics 91

Letters. — 1985. — 12. — Vol. 122, no. 3. — Pp. 264–270. — URL: http: //linkinghub.elsevier.com/retrieve/pii/0009261485805741. [64] Ziegler J., Biersack J., Littmark U. No Title // Stopp. Range Ions Solids. — 1985. — Vol. 1. — Pp. 109–140. [65] Brenner Donald W. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films // Physical Review B. — 1990. — 11. — Vol. 42, no. 15. — Pp. 9458–9471. — URL: http://link.aps.org/doi/10.1103/PhysRevB.42.9458. [66] Tersoff J. New empirical approach for the structure and energy of covalent systems // Physical Review B. — 1988. — 4. — Vol. 37, no. 12. — Pp. 6991– 7000. — URL: http://link.aps.org/doi/10.1103/PhysRevB.37.6991. [67] Wang Binjun, Urbassek Herbert M. Molecular dynamics study of the {\alpha}–{\gamma} phase transition in Fe induced by shear deformation // Acta Materialia. — 2013. — 9. — Vol. 61, no. 16. — Pp. 5979–5987. — URL: http://linkinghub.elsevier.com/retrieve/ pii/S1359645413004710. 92

Рецензии:

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

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

Отзывы:

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