Размещено на http://www.allbest.ru/

Министерство высшего и среднего специального образования республики Узбекистан

Национальный Университет Узбекистана имени Мирзо Улугбека

химический факультет

Кафедра физической и коллоидной химии

Выпускная работа

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

ИБАДУЛЛАЕВА ?ИЛОЛА АЛИШЕР ?ИЗИ

Научный руководители: д.х.н, проф. Акбаров Х.И

ф-м.ф.н. Мухторов А.П

Ташкент-2012

Ўзбекистон республикаси олий ва ўрта махсус таълим вазирлиги

мирзо улу?бек номидаги ўзбекистон миллий университети

Кимё факультети

Физикавий ва коллоид кимё кафедраси

БИТИРУВ ИШИ

КИМЁВИЙ МОДЕЛЛАШДА МОНТЕ КАРЛО УСУЛИ

ИБАДУЛЛАЕВА ?ИЛОЛА АЛИШЕР ?ИЗИ

Илмий ра?барлар: проф. Акбаров Х.И

ф-м.ф.н. Мухторов А.П

Тошкент-2012

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

ГЛАВА I. КВАНТОВОХИМИЧЕСКИЕ ИССЛЕДОВАНИЯ НАНОРАЗМЕРНЫХ КЛАСТЕРОВ КРЕМНИЯ

1.1 Методы исследований электронной и пространственной структуры наноразмерных кластеров

1.1.1 Квантовохимические методы исследований

1.1.2 Методы сильной связи

1.1.2.1 Обзор методов сильной связи

1.1.2.2 Нестандартный метод сильной связи

1.2 Метод Монте Карло в моделировании

1.2.1 Общие теоретические положения метода Монте Карло

1.2.2 Алгоритм Метрополиса

ГЛАВА II. Обзор исследований вакансии в Нанокремнии

ГЛАВА.III Исследование электронной и пространственной структуры вакансии в нанокремнии

3.1 Исследование устойчивости и зарядовых состояний вакансии в димеризованном кластере Si29

3.2 Влияние вакансии на структурные параметры и энергетические характеристики кластеров кремния Si38 и Si38H30

Заключение и выводы

Список цитированной литературы

ГЛАВА I. КВАНТОВОХИМИЧЕСКИЕ ИССЛЕДОВАНИЯ НАНОРАЗМЕРНЫХ КЛАСТЕРОВ КРЕМНИЯ

1.1 Методы исследований электронной и пространственной структуры наноразмерных кластеров

1.1.1 Квантовохимические методы исследования

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

Методы расчета электронной структуры многоатомных систем подразделяются на три уровня сложности: неэмпирические (ab-initio) [1-2], эмпирические и полуэмпирические. Хотя эмпирические модели (межатомные потенциалы) являются с вычислительной точки зрения наиболее эффективными, они не содержат явной информации об электронной структуре, что сильно ограничивает их применимость, например, к неизвестным структурам или структурам с различным типом химической связи.

Неэмпирические расчеты традиционно разделяются на квантовохимические и квантовомеханические подходы. В основе квантовохимических методов лежит широко известный метод Хартри-Фока. Метод Хартри Фока, в котором учитывается взаимодействие электрона с полем, получаемым путем усреднения по всем положениям остальных электронов, основан на решении системы уравнений Хартри (по одному уравнению для каждого занятого одноэлектронного уровня i(r))[1,3]

. (1.1)

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

. (1.2)

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

Степень строгости расчетов, а точнее их сложности, может быть различной: в вычисления с разной степенью полноты можно включить возбужденные электронные конфигурации, а можно ограничиться однодетерминантным приближением Хартри-Фока. Число базисных функций, используемых при расчетах, также может быть различным. Все определяется целью расчета и имеющимися компьютерными ресурсами. Ясно, что чем более высокого уровня расчет, тем более точные результаты могут быть получены с его помощью. В случаях, когда исследуются ряды соединений и важны лишь относительные, а не абсолютные значения энергии и других характеристик, можно не вычислять интегралы, возникающие в схеме расчета, а оценивать их значения на основании экспериментальной информации. Однако использование неэмпирических методов для расчета больших многоэлектронных и многоядерных систем связано с огромными вычислительными трудностями, которые не удается преодолеть даже с помощью современных ЭВМ.

Одним из наиболее широко принятых квантовомеханических методов расчета является метод функционала плотности. Метод основан на предположении, что энергия основного состояния системы взаимодействующих частиц есть однозначный функционал распределения плотности этих частиц (r), и, следовательно, (r) неявно определяет все свойства основного состояния системы. Истинное распределение плотности дается минимумом функционала энергии [2]

Eel((r))/ (r)=0. (1.3)

Решение такого уравнения для плотности (r) было предложено искать в виде плотности невзаимодействующих электронов, находящихся во внешнем самосогласованном поле

. (1.4)

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

. (1.5)

Здесь - эффективный потенциал, который имеет вид

, (1.6)

а обменно-корреляционный потенциал Vxcравен

, (1.7)

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

, (1.8)

где функция определяется свойствами однородного взаимодействующего электронного газа, то задача вычисления плотности электронов в кристалле и его энергии в основном состоянии полностью определяется описанной выше системой уравнений [3].

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

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

Дальнейшие приближения, это - относительно простые и достаточно надежные полуэмпирические методы, в том числе методы сильной связи (МСС) [6-22], в настоящее время стали популярными и удобными методами расчета полной энергии, а также компьютерного моделирования сложных систем. Эти методы позволяют сильно сокращать вычислительное время путем введения подгоночных параметров вместо расчета огромного количества сложных интегралов, одновременно сохраняя квантовомеханический характер расчетов.

В стандартных МСС [4-5] полная энергия системы записывается в виде

, (1.9)

наноразмерный кластер монте карло

где i и Ni - энергия i-го электронного состояния и число электронов в нем. Первый член в (1.9) соответствует энергии отталкивания ядер (или остовов атомов) минус энергия отталкивания электронов; последняя входит дважды во второй член, а последний - энергии изолированных атомов, используемых нередко в качестве подгоночного параметра.

Различие стандартных МСС друг от друга заключается в основном в выборе формул, аппроксимирующих матричные элементы гамильтониана и энергию отталкивания в зависимости от расстояния между атомами. В ранних вариантах МСС [4,6] для этих зависимостей выбиралось простое r-2 масштабирование, которое, будучи удовлетворительным для систем с алмазоподобной структурой, оказалось непригодным для плотноупакованных систем. В настоящее время существуют в основном три улучшенных варианта стандартного МСС, широко используемых для систем преимущественно с ковалентной связью.

Гудвин с сотр. [5] (GSP) показали возможность улучшения переносимости МСС на другие структуры при выборе следующей функциональной зависимости для матричных элементов

. (1.10)

Потенциал отталкивания представлялся в виде суммы по парным потенциалам (rij) с похожей функциональной формой. В данном варианте МСС учитывается взаимодействие только первых ближайших соседей кремния.

Модель Kwon с сотр. [6] (KWON) использует другие значения nc и rc для каждого типа матричного элемента типа s-s, s-p, p-p, различая также - и -взаимодействия. Потенциал отталкивания принимает форму, отвечающую модели погруженного в среду атома,

. (1.11)

В этом варианте МСС радиус взаимодействия обрезается на расстоянии 4.16 Е, поэтому он включает, как правило, взаимодействие с атомами более одной ближайшей координационной сферы.

Третий вариант стандартного МСС (KOH) имеет другую формулу для матричных элементов гамильтониана:

. (1.12)

Потенциал отталкивания имеет похожую функциональную форму с дополнительным фактором Aij, который зависит от локального окружения

, (1.13)

где Zi - эффективное координационное число i-го атома, определяемое в виде:

. (1.14)

Данный вариант МСС в случае алмазоподобной структуры учитывает взаимодействие с атомами вплоть до третьей координационной сферы.

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

Сравнительный анализ указанных вариантов МСС [7] для исследования дефектов в кристаллическом кремнии показал, что увеличение сложности функциональных форм не обязательно приводит к улучшению получаемых результатов. Для большинства статических и динамических свойств ни одна из этих форм не оказалась лучше остальных. Более того, самосогласованный вариант стандартного МСС [5-7] также не приводит кулучшению результатов, например, в случае исследования кластеров кремния.

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

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

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

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

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

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

1.1.2 Метод сильной связи

Для расчета энергии многоатомных систем имеются приближения трех уровней сложности: эмпирические [1-3], полуэмпирические [3-5] и первопринципные (ab-initio) приближения. Хотя эмпирические модели (межатомные потенциалы) являются с вычислительной точки зрения наиболее эффективными, они не содержат явной информации об электронной структуре, что сильно ограничивает их применимость, например, к неизвестным структурам или структурам с различным типом химической связи. В последнее время в связи с повышением мощности компьютеров и разработки новых алгоритмов (таких, как методы с N-скейлингом [5-7], где N - число частиц в системе) расчеты электронной структуры становятся в вычислительном отношении существенно более эффективными. Однако вычислительные затраты методов из первых принципов являются еще достаточно высокими, так что эти методы являются непрактичными для моделирования большинства систем, в особенности их временной эволюции.

Существующие варианты МСС можно разделить на две группы по способу представления полной энергии системы. В методах первой группы [4-6,9-20], которые назовем стандартными МСС, полная энергия системы записывается в виде

, (1.14)

где i и Ni - энергия i-го электронного состояния и число электронов в нем, первый член в (1.14) соответствует энергии отталкивания (энергия отталкивания ядер (или остовов атомов) минус энергия отталкивания электронов; последняя входит дважды во второй член), а последний - энергии изолированных атомов, используемых нередко в качестве подгоночного параметра.

1.1.2.1 Обзор методов сильной связи

Различие стандартных МСС друг от друга заключается в основном в выборе формул, аппроксимирующих матричные элементы гамильтониана и энергию отталкивания в зависимости от расстояния между атомами. В ранних вариантах МСС [9] для этих зависимостей выбиралось простое r-2 масштабирование, которое, будучи удовлетворительным для систем с алмазоподобной структурой, оказалось непригодным для плотноупакованных систем. В настоящее время существуют в основном три улучшенных варианта стандартного МСС, широко используемых для систем преимущественно с ковалентной связью.

Гудвин (GSP) показали возможность улучшения переносимости МСС на другие структуры при выборе следующей функциональной зависимости для матричных элементов

. (1.15)

Потенциал отталкивания представлялся в виде суммы по парным потенциалам (rij) с похожей функциональной формой. В данном варианте МСС учитывается взаимодействие только первых ближайших соседей, за исключением структур типа -олова, где учитывается также взаимодействие вторых соседей.

Модель Kwon с сотр. [5-7] (KWON) использует другие значения nc и rc для каждого типа матричного элемента типа s-s, s-p, p-p, различая также - и -взаимодействия. Потенциал отталкивания принимает форму, отвечающую модели погруженного в среду атома,

. (1.16)

В данном варианте МСС радиус взаимодействия обрезается на расстоянии 4.16 Е, поэтому он включает, как правило, взаимодействие с атомами более одной ближайшей координационной сферы.

Третий вариант стандартного МСС [9] (KOH) имеет другую формулу для матричных элементов гамильтониана:

. (1.17)

Потенциал отталкивания имеет похожую функциональную форму с дополнительным фактором Aij, который зависит от локального окружения

, (1.18)

где Zi - эффективное координационное число i-го атома, определяемое в виде:

. (1.19)

Данный вариант МСС в случае алмазоподобной структуры учитывает взаимодействие с атомами вплоть до третьей координационной сферы.

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

Сравнительный анализуказанных вариантов МСС для исследования дефектов в кристаллическом кремнии показал, что увеличение сложности функциональных форм не обязательно приводит к улучшению получаемых результатов. Для большинства статических и динамических свойств ни одна из этих форм не оказалась лучше остальных. Более того, самосогласованный вариант стандартного МСС [9] также не приводит к улучшению результатов, например, в случае исследования кластеров кремния.

1.1.2.2 Нестандартный метод сильной связи (НМСС)

Предложенный в работе [4-5] вариант МСС, который назовем нестандартным МСС (НМСС), в отличие от стандартных МСС основан на другом выражении для функционала полной энергии и использует относительно небольшое количество параметров. Причем эти параметры имеют определенный физический смысл и связаны с показателями экспоненты слейтеровских атомных орбиталей (АО), которые характеризуют меру протяженности волновых функций электронов в атомах. Кроме того, радиус взаимодействия явно не ограничивается, однако выбранная экспоненциальная зависимость для матричных элементов (см. далее) приводит к их быстрому спаду и исчезновению уже для третьих или четвертых соседей.

В НМСС полная энергия относительно энергии изолированных атомов () 21 представляется в следующем виде:

, (1.20)

где м и н обозначают ядра, а i и j - атомные орбитали (АО); и - экранированные заряды ядер в нейтральных атомах и неточечные ионные заряды; Rмн - межатомные расстояния; H и P - матрицы гамильтониана и электронной плотности (порядков связей) системы; и - полные энергии индивидуальных атомов (ионов), отвечающие соответственно их изолированным состояниям и состояниям во взаимодействующей системе.

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

Электронные энергетические уровни kисследуемой системы, а также коэффициенты разложения волновых функций электронов по АО определяются с помощью решения секулярного уравнения

(1.21)

в предположении об ортогональности АО, как и в отмеченных выше МСС, но самосогласованным расчетом распределения электронов по АО. Отметим, что во всех МСС учитываются только валентные электроны, остальные электроны вместе с ядром составляют остов атома.

Формулы для , , () и выражаются через Rмн, число электронов и на атомных орбиталях (АО), показатели слейтеровских экспонент АО, а также содержат небольшое количество подгоночных параметров

, (1.22)

, , (1.23)

, (1.24)

, (1.25)

=, (1.26)

=, (1.27)

, , (1.28)

где fij - функции взаимной ориентации АО i и j [3,4], и - наиболее вероятное и среднее расстояния между электроном на ядром, n и - главное квантовое число и показатель слейтеровской АО.

Самосогласованные расчеты электронной структуры основаны на пересчете энергии АО в зависимости от степени их заселенности электронами. Для этой цели для валентных электронов используется формула Слетера-Зенера [40] в несколько модифицированном виде

, (1.29)

, , (1.30)

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

Оказалось [8], что формула

(1.31)

является хорошим приближением для чисел экранирования Sij, если для принимается линейная зависимость от зарядового состояния атома, , в виде

=, если ,

=, если , .

При этом значения Sij лежат в достаточно узком интервале значений для всех элементов с s и p валентными электронами и использование даже единственной константы экранирования позволяет достаточно хорошо описывать [21] электронное сродство и до четырех потенциалов ионизации этих элементов.

В [21] последний член (1.20) рассчитывается не с помощью (1.29), (1.30) используется только в расчете (1.22), а следующей простой формулой, выраженной непосредственно через первый потенциал ионизации Iи электронное сродство A атома

. (1.31)

Отметим, что q, как и , измеряется в единицах заряда электрона, т.е. безразмерная величина.

Данный вариант НМСС использует следующие значения параметров: (параметр, используемый в (1.25), а не в (1.31), , и для кремния. Первые три параметра были определены в [21] из условия воспроизведения энергии и длины связи, а также частоты колебаний нейтральной двухатомной молекулы Si2. НМСС даже в такой простой параметризации, включающей три параметра, позволяет объяснить ряд таких нетривиальных эффектов [41-42], как U-отрицательные свойства вакансии в кремнии, термическая миграция собственного междоузельного атома в кремнии и т.д.

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

В НМСС функционал полной энергии представляется в следующем виде:

, (1.32)

где - межядерное расстояние,

, (1.33)

(1.34)

соответственно экранированный ядерный и неточечный ионный заряды; - заряд -го ядра (или ядра вместе с остовными электронами); - наиболее вероятное расстояние между электроном и ядром, nи - главное квантовое число и слейтеровские экспоненты i-ой AO, центрированной на -ом ядре; и - полные энергии индивидуальных атомов в невзаимодействующей и взаимодействующей системах, характеризующихся числами заполнения {} и {} и энергиями {} и {} валентных AO соответственно; и a - подгоночные параметры. Здесь и далее всегда подразумеваются атомные единицы величин.

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

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

AO предполагаются ортогональными и матричное уравнение вида

(1.35)

решается самосогласованно для определения энергетического спектра {} и коэффициентов разложения Сj молекулярных орбиталей (МО) системы по AO. Самосогласованные расчеты основаны на итеративном перерасчете диагональных матричных элементов гамильтониана с использованием зависимости матрицы порядков связей

(1.36)

и заселенности АО от Сj. Здесь k нумерует MO, Nk - заселенности MO, в традиционных МСС равные 2 (за исключением высших занятых MO систем с нечетным числом электронов). Однако НМСС предполагает описание заряженных и возбужденных систем, где одна и более MO может иметь заселенность, меньшую 2.

Диагональные и недиагональные матричные элементы НМСС имеют соответственно вид

(1.37)

и

, (1.38)

где

, (1.39)

, (1.40)

- среднее расстояние между электроном и соответствующим ядром, - угловые функции, протабулированные Слейтером и Костером [6]. В (1.38) знак плюс берется для sp и pp- матричных элементов, а знак минус - для ss и pp- матричных элементов.

Здесь и ниже в формулах , , , a, b, d, e - подгоночные параметры, расстояния измеряются в единицах боровского радиуса, энергии - в хартри, заряды - в единицах заряда электрона.

и определяются следующими формулами [39]

, (1.41)

, (1.42)

где

, , (1.43)

, , (1.44)

, , (1.45)

, (1.46)

, (1.47)

, . (1.48)

Модификация [41] формулы Слейтера-Зернера состоит в использовании произведения квантового числа и эффективного квантового числа в (1.41) вместо квадрата последнего, а также предположения о линейной зависимости эффективного квантового числа от зарядового состояния атомов, так что в пределе водородоподобного иона получаются квантовые числа для атома водорода. Постоянные экранирования (1.44) позволяют с разумной точностью рассчитывать энергии AO для нескольких зарядовых состояний элементов, а также соответствующие потенциалы ионизации и электронные сродства. Для точного воспроизведения нескольких важных энергетических переходов в атомах и ионах вводятся подгоночные параметры через (1.48). Первые два подгоночных параметра в каждой из формул (1.46),(1.47) корректируют универсальные постоянные экранирования (1.44). Эти данные, имеющиеся почти для всех элементов периодической системы, являются базовыми данными НМСС.

Последний член (1.32) рассчитывается с помощью простой формулы, выраженной непосредственно через первый потенциал ионизации Iи электронное сродство A атома

. (1.49)

Отметим, что q, как и , измеряется в единицах заряда электрона, т.е. безразмерная величина.

Найденные из (1.36) числа заполнения АО используются для расчета (формула (1.41)) и неточечных зарядов Q (формула (1.34)), следовательно, для пересчета диагональных матричных элементов (1.37). Итеративное решение уравнения (1.35) прекращается, когда максимальное изменение в двух последовательных итерациях становится меньше наперед заданного малого числа порядка 10-6-10-9. Причем высокая точность используется на этапе параметризации НМСС или при расчете частот колебаний атомов. Во избежание лишних вычислений полная энергия (1.32) рассчитывается только после достижения самосогласования.

Для ускорения сходимости самосогласования в программе расчетов мы реализовали процедуру динамического демпфирования зарядов (ПДЗ). ПДЗ служит для подавления осцилляций зарядов и предотвращения нефизического пераспределения зарядов в системе. Это достигается путем пересчета матрицы порядков связи для текущей итерации Pn, рассчитанной по формуле (16), с помощью следующего алгоритма

, (1.50)

где

,

, .

Поскольку матрицы из двух предыдущих итераций, Pn-1 и Pn-2, будут доступными только при третьей итерации, в первых двух итерациях мы используем постоянное значение параметра демпфирования (0 <a< 1). Далее Pn-1 и Pn-2 заменяются на P*n-1 и P*n-2, начиная с n = 3 и n = 4 соответственно[9].

1.2 метод Монте Карло В моделированиИ

1.2.1 Общие теоретические положения метода Монте Карло

Метод Монте Карло, применяемый в статистической физике, является частным случаем общего метода статистического моделирования, который используют для решения широкого круга задач в различных областях науки. В рамках метода Гиббса термодинамические характеристики вещества получают в результате усреднения по ансамблю, т.е. по совокупности очень большого числа идентичных по природе систем, находящихся в одинаковых внешних условиях и различающихся только по микросостоянию [10-11]. Принято выделять три типа ансамблей, состояние которых задается три типами функций распределения.

В микроканоническом ансамбле (N,V,E) рассматривают замкнутые изолированные системы, в которых фиксированы число частиц N, объем V и полная энергия E. На микроскопическом уровне существует бесконечное число различных способов, или конфигураций, в которых может быть реализовано данное макросостояние. Они могут обмениваться энергией и частицами с окружающей средой. При этом полагают, что рассматриваемая система мала по сравнению с окружающей её системой, и любое изменение характеристик малой системы не сказывается на состоянии большой. Большая система действует как тепловой резервуар или тепловая баня с заданной абсолютной температурой T. В большом каноническом ансамбле(T,V,µ) системы способны обмениваться и энергией, и частицами. Его состояние задается температурой T, объёмом V и химическим потенциалом µ. Расчет термодинамических характеристик, как правило, проводится в рамках канонического ансамбля(N,V,T; N,P,T). Каноническое распределение Гиббса статистическое распределение для систем, содержащих заданное число частиц N, объём V(или давление (P) и способных обмениваться энергией с окружением[12].

Вероятность нахождения системы вмикросостоянии i c энергией Еi рассчитывают по формуле:

wi = (1/Z) exp (-Ei/kT) (1.1)

где k - постоянная Боьцмана, Z- статистическая сумма, сумма по состояниям системы:

Z=У exp (-Ei/kT) (1.2)

В классической статистической термодинамике микросостояние определяется заданием обобщенных координат q и обобщенных импульсов p. Для канонического ансамбля N частиц вероятность иметь значения импульсов в интервале(p, p+? p) значения координат в интервале (q, q+? q) определяется как:

dW(p,q)=(1/Z) exp [-H(p,q)/kT]dpdq (1.3)

Для классической системы статистическая сумма (1.2) заменяется статистическим интегралом

Z= dpdq (1.4)

fi-число степеней свободы молекул сорта i.

Здесь

H(p,q)=+ + UN(q1…..qN) (1.5)

Где m -- масса молекулы. Функцию Гамильтона, отсчитываемую от нулевой энергии молекул, можно представить как сумму энергии внутренних молекулярных движений (электронные состояния, колебания, вращения и т.д.) Нвнутр, энергии поступательного движения центров масс и потенциальной энергии межмолекулярных взаимодействий[13].

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

Q = ? Qпост ? Qвнутр(1.6)

где E0 - энергия молекулы в самом низком энергетическом состоянии, Qпост - статистическая сумма, связанная с поступательным движением молекулы, Qвнутр - статистическая сумма, связанная с внутренними молекулярными движениями.

Подставляя выражение для Гамильтониана в (1.3) и, учитывая условие нормировки, получим выражение для статистического интеграла в виде

Z = (1.7)

где (1.8)

- конфигурационный интеграл. Если статистический интеграл известен, то для рассматриваемой системы можно найти все термодинамические величины; так свободная энергия Гельмгольца системы[14] равна

F = - kT ln Z (1.8)

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

p=S== (1.9)

Конфигурационный интеграл Zконф как функция температуры T, объема V и числа частиц N дает полную статистико-механическую информацию о системе и позволяет по формулам (1.7)-(1.9) рассчитать термодинамические свойства системы.

Невозможность точного вычисления конфигурационного интеграла для реальных систем приводит к необходимости применения новых методов, в которых избегают непосредственного вычисления Zконф. Одним из таких методов расчета является метод Монте-Карло[13-15].

1.2.2 Алгоритм Метрополиса

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

wi = (1.10)

Если в некоторый фиксированный объем помещать случайным образом молекулы, энергия взаимодействия между которыми задается набором потенциальных функций, то в зависимости от конфигурации системы больцмановский множитель exp (- Ui /кТ) может принимать различные значения[11-17]. Некоторые конфигурации дают значительный вклад в канонические средние, а некоторые - практически нулевой (например, когда две частицы сближены настолько, что между ними имеется сильное отталкивание).

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

Среднее по ансамблю от любой физической величины M рассчитывают по формуле

<M>= iexp(-)/= (1/11)

где i - номер конфигурации (среднее берется по всем рассмотренным конфигурациям системы)[16]. Поскольку усреднение (1.11.) проводят по конечному числу конфигураций m со смещенной выборкой, для исключения влияния смещения каждую конфигурацию необходимо брать с весом 1/рi :

<M>? exp(-)/ (3)

Метрополис с соавторами [15-16] предложил в качестве рi взять распределение Больцмана

Рi = exp()/ (1.13)

В результате среднее значение любой физической величины M можно записать в виде

<M>? (1/m) (1.14)

Ансамбль, состоящий из m конфигураций, получают путем задания вероятностей перехода от одной конфигурации к другой. Вероятность перехода от i-й конфигурации к j-й pij считают зависящей от энергий этих конфигураций, а точнее от величины (Uj -Ui)/kT:

Рij = Pji exp[-(Uj -Ui)/kT] (1.15)

Таким образом, строят простые цепи Маркова, т.е. последовательности случайных событий, в которых вероятность определенного события зависит от исхода предыдущего испытания. В соответствии с условием микроскопической обратимости вероятности pij должны удовлетворять условиям:

Рij=exp() Pji=exp() и (1.16)

Обычно полагают, что при

Uj ? Ui Pji = wi (1.17)

и при Uj >Ui Рij = wi exp[-(Uj -Ui)/kT]

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

На практике алгоритм Метрополиса реализуют следующим образом. Пусть заданы потенциал взаимодействия, конфигурация системы (начальное расположение частиц в элементарной ячейке моделирования) и температура Т. Рассчитывают потенциальную энергию системы Ui и вносят случайное изменение в конфигурацию (случайным образом выбирают к -ю частицу в ячейке и смещают ее). При этом энергия системы становится равной Uj. Если Uj < Ui , то считают, что система перешла в новое состояние. Если Uj < Ui, то сравнивают величину exp() со случайным числом оє(0,1). Если о? exp[-(Uj -Ui)/kT], то считают, что система перешла в j-е состояние. Если же о> exp[-(Uj -Ui)/kТ ], то переход в новое состояние не происходит, к -я частица сохраняет свои прежние координаты [18]. При этом j-ую конфигурацию в цепи не учитывают, а рассматривают прежнее расположение частиц, соответствующее энергии U Таким образом, чем больше значение энергии имеет система при случайном изменении конфигурации, тем с меньшей вероятностью она переходит в это состояние.

Очевидно, что расчет средних значений сопряжен с большим количеством вычислений, выполнение которых стало возможным только при появлении быстродействующих ЭВМ[19].

Одно из преимуществ метода Монте-Карло состоит в том, что алгоритм легко адаптировать к любому статистическому ансамблю. Например, при моделировании системы в NPT-ансамбле необходимо периодически изменять объем ячейки, а в уравнении (1.17) вместо разности потенциальных энергий использовать разность энтальпий:

ДH = ДU + PДV - kТ ln(1+ДV/V)N, (1.18)

где P - давление, V - объем ячейки, ДV - изменение объема ячейки.

Усредняя термодинамические функции по конфигурациям на равновесном участке цепи Маркова по уравнению (1.14), легко рассчитать конфигурационную энтальпию и мольный объем :

Нконф = <U> + P<V>, Vm = <V>NA/N (1.19)

Для расчета других термодинамических функций - изобарической теплоемкости СР, изотермической сжимаемости kT и коэффициента температурного расширения б, необходимо проводить моделирования системы при различающихся параметрах состояния, а затем находить конечные разности[20]. Оценку указанных термодинамических функций можно сделать и по результатам одного моделирования, используя флуктуационные формулы:

kT==),

==),

Для решения большинства задач вполне достаточно проведения вычислений в каноническом ансамбле. Однако если требуется, то можно использовать и большой канонический ансамбль [21, 22].

Глава II. Обзор исследований вакансии в Нанокремнии

2.1 УСТОЙЧИВЫЕ КОНФИГУРАЦИИ И ЗАРЯДОВЫЕ СОСТОЯНИЯ ОСНОВНЫХ СОБСТВЕННЫХ ДЕФЕКТОВ В КРЕМНИИ

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

Вакансия может существовать в пяти зарядовых состояниях (V++, V+V0, V-. У--) [23]. Вакансия очень подвижна и ее непосредственно удалось идентифицировать только в кремнии р-типа. Она отжигается в температурном интервале 150-180°С с энергией активации "0.33 эВ. Вакансия считается более подвижной в кремнии п-типа (энергия активации -0.18 эВ). Неэмпирические расчеты в рамках ФЛП энергий миграции вакансии дают значения ~1 эВ. Отметим, что ошибки ФЛП-расчетов энергии миграции дефектов составляет 20.5 эВ [24]

Ян-теллеровское искажение решетки, снимая энергетическое вырождение, приводит к расщеплению трехкратного уровня на более низкий однократный и более высокий двукратный уровни. Искажение является тетрагональным (атомы - первые соседи попарно сближаются друг к другу) и имеет группу симметрии D2d.

Влияние ян-теллеровских искажений специально рассматривалось в работах. Поскольку метод ФГ не очень пригоден для адекватного описания искажений решетки вокруг дефектов, этот недостаток старались преодолеть путем введение дополнительных орби талей на расстоянии 0.5 А от каждого из четырех соседей вакансии. Исследовалось положение локализованных уровней вакансии в зависимости от трех типов искажений: а) тетрагонального, б) полносимметричного и в) комбинации первого со вторым. Оказалось, что тетрагональные искажения или их комбинации с полносимметричным в случае смещения атомов во внутрь приводят к расщеплению уровняt2 симметрии на величину порядка полуширины запрещенной зоны кремния.

Одним из наиболее интересных результатов, полученных за последние 10-15 лет в теории дефектов полупроводниковых кристаллов, является предсказание отрицательного значения эффективной энергии взаимодействия электронов (Ч) для вакансии в кремнии (И-отрицательный характер вакансии).

Согласно этим работам три зарядовых состояния вакансии V0, V+ и V++ образуют так называемую "U-отрицательную систему Андерсона" [23-25]. Электронный переход, обусловленный изменением уровня Ферми, происходит между состояниями V0 и V++, тогда как V+ не является стабильным состоянием низкой энергии. На первый взгляд переход второго электрона на однократный занятый уровень должен увеличить энергию системы за счет электрон-электронного отталкивания. Однако наличие второго электрона благоприятствует энергопонижаюшим структурным искажениям. Выигрыш энергии, обусловленный последними, с избытком компенсирует энергию электрон-электронного отталкивания. Таким образом, наблюдается интересное явление: два электрона или одновременно садятся на уровень вакансии, или одновременно покидают ее. Полученные в /5/ значения энергии переходов Е(0/+)=0.05 эВ и Е(+/++)=0.13 эВ находятся в прекрасном согласии с экспериментом [26].

В работе [27] представлены вычисления стабильности и оптических свойств наноматериалов полупроводника с восстановленными поверхностями методами функционала плотности и квантового Монте Карло. Предсказывает относительную стабильность кремниевых нано-структур с реконструированными и нереконструированными поверхностями, и показывает, что своеобразная трансформация поверхностных конфигураций для чрезвычайно кривых поверхностей резко уменьшают оптические промежутки и уменьшают времена жизни экситонов. Эти предсказания обеспечивают объяснение как изменения в спектрах фотолюминесценции синтезированных коллоидным способом нано-частицы, так и наблюдаемые глубокие уровни в запрещенной зоне пористого кремния. На работе демонстрирует, что электрохимическим травлением получены насыщенные водородом кластерыSinHx с n больший чем 20, в пределах семьи дискретных размеров. Эти размеры - 1.0 (Si29), 1.67 (Si123), 2.15, 2.9, и 3.7 нм в диаметре. Они характеризуют частицы через прямое электронное отображение, возбуждения и эмиссионной оптической спектроскопии, и коллоидной кристаллизации. Измерены зоны запрещенных уровней и полосы эмиссии. Наименьшие четыре - ультраяркие синие, зеленые, желтые и красные люминесцентные частицы. Наличие дискретных размеров и отличной эмиссии в красном, зеленом и синем (RGB) диапазонах полезна для биомедицинской маркировки, RGB дисплеев, и флэш-памяти. Посредством крупномасштабных квантовых моделирований мы исследуем формирование и закрепление кластеров вакансии Vn в кремнии для n ? 35. Мы показываем, что существуют различные модели роста и что взаимодействие между аргументами энергии и топологии определяет самые устойчивые агрегатами. Используя метод возбуждения мультистолкновениями (multicollision excitation method) было исследовано, диссоциация ионов кремниевых кластеров, содержащих до 70 атомов. Кластеры с 19-35 атомами разделяются главным образом потерей таких частей, как Si10. Большие кластеры фрагментируются потерейSi6. Для кластеров с n> 60, потеря индивидуального атома Si (как происходит для объемного кремния) начинает становиться важным каналом диссоциации. Эти изменения в наблюдаемых продуктах могут быть приписаны изменениям в энергиях связи кластеров. Когда большие кластеры чрезвычайно возбуждены, они диссоциируют последовательной потерей единиц Si6, пока они не достигают количества атомов ~30. Тогда они распадаются на фрагменты, содержащие 6-11 атомов. Энергии диссоциации были оценены от экспериментальных результатов. Существенные изменения в энергиях диссоциаций наблюдены для кластеров с 25 атомами. Для кластеров с n> 25 энергии диссоциации увеличиваются гладко[27]. Структуры самых низко-лежащих по энергии кластеровSin (n = 20-40 и 45) получены с генетическим алгоритмом и в рамках не-ортогональной модели сильной связи (NTB). Мы находим, что атомы организованы в две оболочки. Внешний формирует искаженную клетку, а внутренний содержит один атом или несколько основных атомов (ядро кластера). Геометрия кластеров Sin изменяются на вытянуто-сферические от n = 24 и изменяются на сплющенно-сферические структуры от n = 29 и изменяются на почти сферические структуры от n = 32. С увеличением размеров кластеров n энергии связи в этом диапазоне увеличиваются гладко, и никакое магическое число не найдено в этом ряду кластеров. Все профили кривых плотностей состояний (DOS) очень подобны. Между профилями соседних кластеров встречаются только небольшие различия. Ширина запрещенной зоны уменьшается с увеличением размера кластеров. Мы вычислили энергии поглощения и эмиссии нанокристаллов кремния в качестве функции размера и пассивантов поверхности, используя как теорию функционала плотности, так и квантовые вычисления Монте Карло. Мы нашли, что ионные перестановки и релаксации электронов, происходящие после поглощения и эмиссии чрезвычайно чувствительны к химии поверхности. В частности нанокластеры с одинаковыми размерами и одинаковыми промежутками поглощения могут показать поразительно различные энергии эмиссии. Наши результаты обеспечивают объединенную интерпретацию нескольких недавних измерений, в которых наблюдали значительно различные энергии эмиссии от кластеров с одинаковыми размерами. Наши вычисления также показывают, что комбинация измерений поглощения и эмиссии может обеспечить мощный инструмент для идентифицирования как размеров, так и пассивантов поверхностей нанокристаллов полупроводника[28]. Метод функционала все-электронной плотности использовалась, чтобы вычислить структуры и энергии кремниевых/водородных кластеров с 148 атомами. Колебательные частоты вычислены для тех кластеров с меньше чем 75 атомами. В дополнение к пассивированным водородом кластерам, основанным на структурах объемного кремния, мы изучаем структуры, включающие вакансии, дивакансии, и дополнительные атомы H. Результаты - по сравнению с более ранней работой и обеспечивают характерные признаки вибрации, которые должны помочь в интерпретации размеров (таких, как инфракрасных спектров) гидрогенизированного кристаллического и аморфного кремния. Предложены и изучены с применением ab inito вычислений эллипсоидальные квантовые точечные структуры. Сравнение энергии формирования со сферическими партнерами показывает маленькие различия. Однако, электронные и оптические свойства явно различны. Поскольку электронные состояния имеют различную чувствительность к изменениям соотношений геометрических размеров, отдельные предельные орбитали отодвигаются в запрещенную зону, и общее перемещение уровней, позволяет учитывать поглощение в видимом диапазоне. Наш результат показывает важность и размера и формы при формировании запрещенной зоны кремниевых квантовых точек[29]. В текущей работе взаимодействие этих 60 степеней перетасовывает дислокацию с группой вакансии под прикладным, стригут напряжение в кремниевом кристалле, изучен через молекулярный метод динамики. Stillinger-Weber (КОРОТКОВОЛНОВЫЙ) потенциальный и зависимый от окружающей среды межатомный потенциал (EDIP) используются, чтобы вычислить межатомные силы. Результаты моделирования показывают, что в низком стригут напряжение, дислокация прикреплена группой вакансии. С уровнем напряжения, увеличенным к определенной критической сигме ценности (l), дислокация может преодолеть скрепление и пройти. Найдено, что сигма (l) достигает своего максимума в температуре перехода, которая является приблизительно 350 K. Также показанный в моделированиях обобщенное разобщение дислокации, что 60 дислокаций степеней, взаимодействуя с группой вакансии, могут привести к 30 степеням и 90 степеням частичные дислокации, когда прикладные стригут пределы напряжения другая критическая сигма ценности (h). Две проистекающих частичных дислокации отделены свойственной ошибкой укладки. В отличие от сигмы (l), сигма (h) продолжает уменьшаться в температурах выше чем 400 K и остается константой в более низких температурах. Электронные структуры кремниевых кластеров с моновакансией и di-вакансией вычислены с использованием дискретного вариационного (DV)-X [альфа] молекулярный орбитальный метод. Результаты показывают, что глубокие уровни энергии произведены в середине энергетического кризиса в сотрудничестве с повисшими облигациями из-за введения вакансии. Мы предлагаем новый механизм устранения глубоких уровней энергии, заменяя вакансиями с атомами или парой атомов для заявлений в Сицзяне ультрабольшие чешуйчатые интегральные схемы. Один из лучших кандидатов как замена вакансии - пара объединенных промежуточных атомов бора[30]