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

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НАНОРАЗМЕРНЫХ СТРУКТУР НА ОСНОВЕ ОКСИДОВ

МЕТАЛЛОВ

Специальность 05.13.18 — математическое моделирование, численные методы и комплексы программ

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

Райк Алексей Владимирович

3 О МАЙ 2013

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

Санкт-Петербург — 2013

005060244

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

доктор физико-математических наук, профессор Егоров Николай Васильевич

доктор физико-математических наук, профессор Андрианов Сергей Николаевич (СПбГУ, факультет ПМ-ПУ)

кандидат физико-математических наук Кримская Ксения Александровна (ОАО «Ленэнерго->)

Московский физико-технический институт (государственный университет) (МФТИ, Москва)

00

Защита состоится «19» июня 2013 г. в 14 часов на заседании совета Д 212.232.50 по защите диссертаций на соискание ученой степени кандидата наук, на соискание ученой степени доктора наук при Санкт--11е гербургском государственном университете по адресу: 198504, г. Санкт-Петербург, Петродворец, ул. Ульяновская, д. 1. НИИ Физики им. В. А. Фока, ауд. 116

С диссертацией можно ознакомиться в Научной библиотеке им. М. Горького Санкт-Петербургского государственного университета по адресу: 199034, Санкт-Петербург, В.О., Университетская наб., 7/9. Автореферат размещен на сайте www. spbu. ru

Автореферат разослан «lfc»_-2013 г.

Ученый секретарь диссертационного совета доктор физ.-мат. наук, профессор

Научный руководитель:

Официальные оппоненты:

Ведущая организация:

Г. И. Курбатова.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

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

1. Разработана методика моделирования кластерных фрагментов поверхности.

2. Рассчитаны энергетические характеристики поверхностных слоев.

3. Исследовано взаимодействие молекулы НгО с поверхностью кристаллов MgO и ZnO.

4. Найдены потенциалы взаимодействия молекулы воды с поверхностями кристаллов.

Методами исследования являются методы математического и компьютерного моделирования, численного эксперимента, методы квантовой механики и физики твёрдого тела.

Положения, выносимые на защиту:

• Кластерные модели фрагментов поверхности с конечным числом атомов, исследуемые в программном пакете Gaussian.

• Способ расчета электронных характеристик, определяющих специфику поверхности кристаллов MgO и ZnO, а также их адсорбционные и каталитические свойства.

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

• Новая математическая модель механизма адсорбции молекулы воды на кристаллических поверхностях М^О и 2пО.

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

• Комплекс программ, реализующий созданную математическую модель.

Научная новизна. Впервые исследованы кластерные модели поверхностей, содержащие до 100 атомов поверхностного слоя. Найдены энергетические уровни кристаллических структур, определены значения зарядов на атомах, а также полная энергия систем.

Предложены новые структуры полупроводников с разными характеристиками запрещенной зоны на основе гранецентрированной кубической кристаллической решетки М^О и ZnO с изменяемым процентным соотношением ионов магния и цинка.

Исследован механизм адсорбции молекул воды на различных кристаллических поверхностях М§0 и ЪпО. Показано, что при взаимодействии с поверхностью М^О энергетически выгодным является образование водородной связи молекулы воды с атомом кислорода поверхности, тогда как адсорбция воды на поверхности кристалла 2пО обусловлена электростатическим взаимодействием.

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

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

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

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

вых задачах. Результаты расчетов соответствуют приведенным в литературе экспериментальным данным, а но предложенным структурам новых материалов подтверждены результатами моделирования на тестовых задачах.

Апробация работы. Результаты диссертационной работы докладывались и обсуждались на следующих конференциях: 40-я международная конференция студентов и аспирантов «Процессы управления и устойчивость» (СПб, СПбГУ, факультет ПМ-ПУ, 2009 г.); 41-я международная конференция студентов и аспирантов «Процессы управления и устойчивость» (СПб, СПбГУ, факультет ПМ-ПУ, 2010 г.); 43-я международная конференция студентов и аспирантов «Процессы управления и устойчивость» (СПб, СПбГУ, факультет ПМ-ПУ, 2012 г.); всероссийская конференция «Устойчивость и процессы управления» (КСР) им. В. И. Зубова (СПб, СПбГУ, 2010 г.); а также на научных семинарах кафедры моделирования электромеханических и компьютерных систем факультета прикладной математики - процессов управления СПбГУ.

Публикации. Основные результаты диссертационной работы отражены в 6 публикациях, список которых приведен в конце автореферата, в том числе 2 работы опубликованы в изданиях, рекомендованных ВАК Министерства образования и науки РФ.

Структура и объем работы. Диссертационная работа изложена на 133 страницах машинописного текста и состоит из введения, четырех глав, заключения и библиографического списка литературы, включающего 139 наименований. Работа содержит 31 рисунок, 19 таблиц и 2 приложения на 13 листах.

СОДЕРЖАНИЕ РАБОТЫ

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

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

Математическая модель основана на нерелятивистском стационарном уравнении Шредингера для многоэлектронной волновой функции Ф:

ЯФ(гь г2, ..., гдг, Нь Кг, • • • > Н-м) = ЕФ(п, г2, ..., тц, Кь Лг, ..., Км),

где Е — полная энергия системы, Н — оператор Гамильтона для системы, состоящей из N электронов и М ядер, Ф(г1, ..., г^, Иг, И-г, • • •, П-м) — молекулярная волновая функция. При отсутствии внешних магнитных и электрических полей Гамильтониан Н в системе атомных единиц имеет следующий вид:

N 1 М 1

^ 2 1 2 Ма

¿=1 а=1

N М _ N N . ММ

-ЕЕ^+ЕЕ^ + ЕЕ

ZaZp

i=lj>irii а=1 (3>с Ra0

Здесь Ма — масса ядер, Za — заряд ядра, Ra — координаты ядер, г* — координаты электронов; греческие индексы относятся к ядрам, а латинские к электронам; Г и V являются операторами кинетической и потенциальной энергии соответственно.

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

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

Ftfii =£i<Pi,

где F — оператор Фока, £, — собственные значения. Метод решения уравнений Хартри—Фока называется методом самосогласованного поля, поскольку предполагается, что каждый электрон движется в поле, создаваемом остальными электронами.

Альтернативой является подход, основанный на использовании теории функционала электронной плотности (Density Functional Theory, DFT). Его идея заключается в том, что распределение электронной плотности р определяет свойства атомов и молекул:

+оо +оо

p(r)=N J ■■■ J |Ф(хь Х2, ...,XN)\2 dX2..-dXN, (1)

—оо —ос

Хм — пространственно-спиновые координаты N-го электрона, р(г) — плотность вероятности нахождения любого из N электронов с произвольным спином в состоянии, описываемом волновой функцией Ф. Знак интегрирования в (1) также подразумевает суммирование по спинам всех электронов. Электронная плотность является неотрицательной функцией трех про-

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

г—>оо ц

N.

V

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

Согласно теореме Хоэнберга—Кона, энергия системы Ео в основном электронном состоянии однозначно определяется электронной плотностью основного состояния ро■ Рассмотрим систему невзаимодействующих электронов, движущихся в некотором статическом потенциале г^(г), который включает в себя внешнее поле г>сх1 (г) и электрон-ядерное притяжение

= (2)

Представим функционал энергии в виде:

Е(р) = ¡ер{т) и(г)йУ + и I + Ехс(р).

V к У1

Первый член в правой части равенства описывает притяжение электронов к ядрам и взаимодействие с внешним полем, второй — энергию кулонов-ского межэлектронного взаимодействия (г* характеризует положение г-го электрона, г,- — электрона,

с!Уг — элемент объема конфигурационного пространства всех электронов системы за исключением г-го электрона, ¿У^ — за исключением j-ro электрона), соответствующую кулоновскому потенциалу

(3)

VI

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

Применение вариационного принципа к функционалу Е(р) при дополнительном условии

p(r) = 2¿|^(r)|2, ¿=i

где <fii(r) — орбитали Кона—Шэма, приводит к системе одноэлектронных интегро-дифференциальных уравнений для орбиталей, называемых уравнениями Кона—Шэма, которые можно записать следующим образом:

¥>n(r¿) =£„¥>п(г»), П = 1,...,ЛГ/2, (4)

где выражения для потенциалов задаются формулами (2) и (3).

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

Проведенный анализ квантово-механических методов позволил предложить оригинальную методику решения поставленных в диссертационной работе задач. На первом этапе использовался метод силового поля для получения оптимизированной геометрии исследуемых кластеров, на этом этапе оптимизация геометрии проводилась методом Ньютона—Рафсона. Под оптимизацией геометрии в работе понимается поиск конфигурации системы, соответствующей минимуму энергии. На втором этапе полученная геометрия использовалась в качестве начальной для квантово-механических расчетов. Полная энергия Е систем находилась квантово-химическими методами Хартри—Фока и одной из реализаций метода функционала электронной плотности B3LYP с представлением волновой функции в валентно-расщепленном базисе 6-31G. Энергия связи молекулы Н20 с поверхностью определена как разность полной энергии системы и энергий изолированных поверхностей и молекулы воды:

АЕСВЯЗЛ = -Ё-поверхность+НгО (¿'поверхности "Ь Ец^о) ■

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

Компьютерное моделирование и численные эксперименты были проведены с использованием программных пакетов: Gaussian, Chemcraft, Open Babel, Maple, Origin, Gnuplot.

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

пг

-7T~Vi + + J(r¿) + 2ra

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

В процессе отработки методики построения кластерной модели были проведены расчеты моделей с разным количеством атомов в поверхностном слое и разным количеством слоев. В частности, рассматривались кластерные структуры одного слоя размерностью [3 х 3]i, [5 х 5]i, [7 х 7]ь двойного слоя [3 х 3]2, [5 х 5]2, [7 х 7]2, а также [9 х 9]i, [9 х 9]2, [9 х 9]4. Минимальное количество атомов поверхностного слоя равнялось 9, максимальное — 324.

Для оптимизации геометрии на первом этапе использовался метод силового поля. В результате для фрагментов поверхности MgO были получены кластеры с параметром решетки ~ 2 А. Плоская структура оптимизированных кластеров не сохранялась, наблюдался выход атомов кислорода над поверхностью.

В качестве начальной геометрии для методов Хартри-Фока и B3LYP использовалась структура поверхности, полученная методом силового поля. Независимо от выбранной начальной геометрии результатом оптимизации являлась поверхность с межатомными расстояниями до ближайших соседних атомов порядка 2,0 А - 2,1 А. Данный результат соответствует известному значению параметра кристаллической решетки MgO1.

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

На основе произведенных расчетов для дальнейших исследований была выбрана двухслойная кластерная модель размерностью [7х7]2. В работе показано, что для фрагментов поверхности [3x3] и [5x5] существенное влияние оказывают эффекты на реакционно более активном краю. В центре структуры [7x7] влияние краевых эффектов ослабляется. Расчеты модели [9x9] давали результаты, сравнимые с результатами, полученными для структуры [7x7]. В связи с этим, с целью экономии временных и вычислительных ресурсов, было решено использовать модель с меньшим количеством атомов и слоев.

В рамках кластерной модели было исследовано взаимодействие двух поверхностей кристалла MgO. Решение этой задачи может являться дополнительным критерием адекватности модели. В качестве моделей поверхности был использован оптимизированный геометрически кластер размерности [5 х 5]2. При исходном расстоянии 2,8 А между поверхностями в процессе оптимизации геометрии они начинали сближаться и смещаться друг относительно друга, образуя монокристалл с межатомным расстоянием 2,1 А, что согласуется с известным экспериментальным значением. Этот результат также подтверждает корректность предложенной в диссертационной работе модели.

Проведены расчеты ширины запрещенной зоны кластеров, состоящих из двух слоев каждый размерностью [7х 7]i, для ZnO и MgO, а также кристал-

1Whited R., Flaten С., Walker W. Exciton thermoreflectance of MgO and CaO // Solid State Communications. - 1973. - Vol. 13. - P. 1903-1905.

Рис. 1. Диаграммы энергетических уровней £п (эВ); а - Mg0/Mg01 Ь - 7лЮ/2пО, с - ЪпО/Мф

ла смешанного типа, в котором один слой представлен поверхностью [7 х 7]ь а другой — ЪъО [7 х 7]х- Результаты расчетов приведены на рисунке (1). '

Ширина запрещенной зоны позволяет определять полупроводниковые свойства материала. Для кристалла ЪпО типа вюрцита (гексагональная сингония) эксперимент дает значение ширины запрещенной зоны 3,37 эВ, для кристалла гпО кубической структуры расчет дает значение 2,45 эВ (рис. 1, Ь). Ширина запрещенной зоны в кристалле — 3,72 эВ (рис. 1, а), это подтверждает его свойства как диэлектрического материала. Композиция в кристаллической решетке атомов и Ъп позволяет получить значение 2,57 эВ (рис. 1, с). Эта величина характеризует новый материал, проявляющий полупроводниковые свойства. Варьируя процентное соотношение ионов гп, можно получать полупроводники с различной ширинои запрещенной зоны. Такие кристаллы могут использоваться при конструировании низкопороговых лазеров, работающих при комнатной температуре. Частота излучения зависит от соотношения ионов магния, цинка и кислорода в кристаллической решетке.

В третьей главе приводятся основные результаты моделирования адсорбции молекулы воды на кристаллах МёО и гпО. При моделировании фрагментов поверхности использовались кластерные модели, описанные в

предыдущей главе.

Супер молекулярная модель с молекулой Н20 над поверхностью кристалла позволяет учесть взаимодействие всех электронов и ядер и конкури-

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

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

Таблица 1. Межъядерные расстояния Я (А), валентные углы Н1-О-Н2 (градусы), заряды на атомах q (а.е.) и волновые числа V (см-1) для молекулы воды над поверхностью кристалла MgO, рассчитанной по кластерной модели

Параметр И 2 ^изолированная Н20 над MgO [7 х 7Ь

R(Mg-O) 2,371

H(0„-Hx) 1,696

ЩЩ-О) 0,9760 1,020

R(G-H2) 0,9760 0,9743

Hi-O-Hj, 108,304 111,502

4„, 0,359 0,525

Чо -0,718 -0,750

Чн, 0,359 0,415

fl 3781 3738

^эксперимент 3756

3616 3138

Эксперимент 3657 3369

1619 1611

^..сп™-,.. 1595 1650

Моделирование позволило найти прочность связи молекулы воды с поверхностью MgO для всех рассмотренных молекулярных кластеров, значение которой лежит в диапазоне 15-22 ккал/моль. Это хорошо совпадает с эмпирической оценкой энергии связи, выполненной методом корреляций

2Foster М., Furse М., Passno D. An FTIR study of water thin films on magnesium oxide // Surface Science. - 2002. - Vol. 502-503. - P. 102-108.

3Wu M. C., Goodman D. W. Acid/base properties of MgO studied by high resolution electron energy loss spectroscopy // Cat. Lett. - 1992. - Vol. 15. - P. 1-11.

4 Wang Y., Nguyen H. N., Truong T. N. Mechanisms of and effect of coadsorption on water dissociation on an oxygen vacancy of the MgO(lOO) surface // Chem. Eur. J. - 2006 - Vol 12. - P. 5859-5867.

J

; / i \ч

Рис. 2. Молекула воды образует водородную связь с фрагментом поверхности

[7x7] 2

(14-17 ккал/моль). Прочность связи в точке минимума полной энергии для фрагмента поверхности [7 х 7]2 составляет 22 ккал/моль.

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

При рассмотрении адсорбции молекулы воды на грани (100) поверхности кристалла ZnO использовалась предложенная методика расчетов. В качестве модели поверхности выбиралась двухслойная кластерная модель [7 X 7]2, основанная на гранецентрированной кубической решетке. Расстояния Zn-O, полученные в результате предварительной оптимизации геометрии кластера [7 х 7]2, составляют 2,38 А-2,41 А. Процесс адсорбции рассматривался над фиксированной поверхностью с вышеуказанными свойствами. Анализ движения молекулы в поле поверхности свидетельствует о существенном отличии механизма адсорбции молекулы Н20 на поверхности MgO и ZnO. Для кристалла оксида цинка не образовывалась водородная связь с кислородом поверхности, как в случае MgO. То есть поверхность кристалла ZnO является более гидрофобной, это косвенно подтверждается и экспериментальными данными7. Энергия связи молекулы воды с поверхностью кристалла ZnO в точке равновесия составляет 13 ккал/моль.

Результаты расчета для поверхности ZnO представлены в табл. (2).

SScamehorii С A., Harrison N. M., McCarthy M. I. Water chemistry on surface defect sites: Chemidissociation versus physisorption on MgO(OOl) // J. Chem. Phys. - 1994. - Vol. 101. - P. 1547-1554.

6Abriou D., Jupille J. Self-inhibition of water dissociation on magnesium oxide surfaces //

Surface Science. - 1999. - Vol. 430. - P. 527-532.

7Tang L., Zhou В., Tian Y. et al. Synthesis and surface hydrophobic functionalization ot ZnO nanocrystals via a facile one-step solution method // Chem. Engin. J. — 2008. — Vol. 139. — P. 642-648.

Таблица 2. Межъядерные расстояния Л (А), валентные углы Н1-О-Н2 (градусы) и заряды на атомах q (а.е.) молекулы воды над поверхностью кристаллов 7пО, рассчитанной по кластерной модели фрагмента поверхности

Параметр Н20из0лнр0ванная Н20 над гпО [7 х 7] 2

И(Н1-0) 0,9759 0,9812

ЩО-Н2) 0,9759 0,9714

Н1-О-Н2 108,304 108,254

Ч„, 0,359 0,445

<Ь -0,718 -0,680

Чн, 0,359 0,439

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

В четвертой главе предложена методика расчета потенциалов взаимодействия молекулы воды с различными поверхностями. Потенциалы вычислялись двумя способами: методом силового поля и методом функционала электронной плотности с гибридным потенциалом ВЗЬУР.

Потенциалы взаимодействия с поверхностью были вычислены методом силового поля для разных ориентаций молекулы воды, приближающейся по нормали к поверхности. В работе показано, что потенциалы существенно зависят от анизотропии системы, определяемой ориентацией молекулы воды над поверхностью. Найдены потенциалы взаимодействия Н2О с поверхностями М^О и ХпО методом ВЗЬУР для наиболее энергетически выгодной конфигурации, когда молекула воды располагалась в плоскости, параллельной плоскости поверхности. Следует отметить, что при этом подходе не рассматривалась возможность образования водородной связи молекулы воды с поверхностью. Вид потенциалов приведен на рисунке (рис. 3). Точка минимума потенциальной кривой соответствует равновесному расстоянию адсорбируемой молекулы до поверхности. В случае взаимодействия с поверхностью 2пО это расстояние больше (кривая, соответствующая ZnO, сдвинута вправо по оси ж), потенциальная кривая для гп0-Н20 имеет более плоский минимум и соответствует меньшему значению энергии связи.

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

Рис. 3. Сравнение потенциалов B3LYP, полученных для Mg0-H20 и Zn0-H20

где R — расстояние между молекулой воды и поверхностью, а А, В, С, D — варьируемые константы.

Для расчета констант А, В, С, D была составлена и решена аналитически система из четырех нелинейных уравнений с четырьмя неизвестными. Данная система имеет не единственное решение. В работе был отобран набор констант, соответствующий их физическому смыслу.

Вид аппроксимирующей функции Бакингема определяется подстановкой R и ip(R) в полученные аналитические выражения для констант. Чтобы учесть физический смысл констант и задать границы их изменения, были выполнены расчеты численными методами в пакетах Gnuplot и Origin. В них реализован нелинейный метод наименьших квадратов Марквардта— Левенберга, являющийся комбинацией методов наискорейшего спуска по градиенту и Гаусса—Ньютона.

Найденные значения констант функции Бакингема, определяющей потенциалы взаимодействия молекулы воды с поверхностями кристаллов MgO и ZnO, представлены в табл. (3):

Таблица 3. Значения констант А, В, С, D систем MgO HaO и ZnO-HaO

Константа Gnuplot Origin

Система Mg0-H20

А -1,76 ± 0,08 -1,75 ± 0,12

В -1,26 ± 0,10 -1,26 ± 0,09

С -0,21 ± 0,04 -0,21 ± 0,02

D 0,81 ± 0,04 0,81 ± 0,02

Система Zn0-H20

А -1,40 ± 0,09 -1,39 ± 0,10

В -2,13 ± 0,06 -2,13 ± 0,07

С -0,23 ± 0,04 -0,23 ± 0,05

D 0,81 ± 0,02 0,80 ± 0,06

Итерационный процесс сходится к одному из минимумов. Большое влияние на конечный результат оказывают начальные значения параметров. Начальной точкой для полученных выше констант служили следующие значения:

А = В = С=-1, D = 1.

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

Проведено исследование скорости расчетов кластерных моделей с адсорбированной молекулой в зависимости от выделенных оперативной и виртуальной памяти, а также от количества ядер процессора и версии программного комплекса Gaussian. Полученные результаты позволяют сделать следующие выводы. Gaussian-09 выполняет расчет для системы Mg0-H20 размерности [5х5]2 существенно быстрее (3 дня 6 часов 37 минут), чем Gaussian-03 (5 дней 0 часов 30 минут). К существенному уменьшению времени расчетов приводит выделение до 4 ядер. Так, масштабируемость при использовании 4 ядер составляет 3,08, при использовании 8 ядер — 4,21. Наиболее ресурсоэффективным является использование 4 ядер для одной задачи.

В заключении сформулированы основные результаты диссертации:

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

2. Представлены результаты расчета геометрии поверхностного слоя и распределения электронной плотности на поверхности кристаллов MgO и ZnO.

3. Предложена новая структура полупроводника на основе изоморфизма кристаллических решеток MgO и ZnO, для которого ширина запрещенной зоны может варьироваться в зависимости от процентного содержания ионов магния и цинка.

4. Установлено принципиальное различие в механизмах адсорбции молекулы Н20 на поверхности MgO и ZnO: молекула воды образует водородную связь с кислородом поверхности MgO вплоть до диссоциации на поверхностных уступах с образованием гидроксильнош иона; взаимодействие с поверхностью ZnO носит электростатический характер.

5 Рассчитаны потенциалы взаимодействия молекулы воды с поверхно-' стями кристаллов М80 и гпО. Полученные значения потенциалов аппроксимированы аналитической функцией Бакингема.

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

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ Публикации в изданиях, рекомендуемых ВАК РФ

1 Райк А. В., Бедрина М. Е. Моделирование процесса адсорбции воды на поверхности кристаллов // Вестн. С.-Петерб. ун-та Сер. 10: Прикладная математика, информатика, процессы управления. 2011 Вып. I. О. ы п.

2 Райк А. В., Егоров Н. В., Бедрина М. Е. Моделирование потенциалов ' межмолекулярного взаимодействия // Вестн. С.-Петерб. С£Р; Ш:

Прикладная математика, информатика, процессы управления. 2012 Вып.

С. 7!>-87.

Публикации в других изданиях

3 Райк А В., Бедрина М. р. Компьютерное моделирование процесса адсорб-' ции воды на поверхности кристалла МйО // Процессы управления и устойчивость: Труды ХЬ международной конференции аспирантов и студентов / под ред. Н. В. Смирнова, Г. Ш. Тамасяна. СПб.: Издат. дом С.-Петерб. гос. ун-та, 2009. С. 247-252.

4 Райк А В. Моделирование потенциалов межмолекулярного взаимодействия ' // Процессы управления и устойчивость: Труды ХЫ международной конференции аспирантов и студентов / под ред. Н. В. Смирнова, Г. Ш. Тамасяна. СПб.: Издат. дом С.-Петерб. гос. ун-та, 2010. С. 205-210.

5 Райк А В. Аппроксимация потенциалов межмолекулярного взаимодействия ' // Процессы управления и устойчивость: Труды ХШ1 международной конференции аспирантов и студентов / под ред. А. С. Ерёмина, Н. В. Смирнова. СПб.: Издат. дом С.-Петерб. гос. ун-та, 2012. С. 183-188.

6 Райк А В., Бедрина М. Е. Моделирование потенциалов взаимодействия адсорбированных молекул с поверхностями кристаллов // Устойчивость и процессы управления. Всероссийская конференция, посвященная 80-летию со дня рождения В. И. Зубова, СПб.: ВВМ, 2010. С. 128-129.

Подписано к печати 07.05.13. Формат 60x84 Лв . Бумага офсетная. Гарнитура Тайме. Печать цифровая. Печ. л. 1,00.

Тираж 100 экз. Заказ 5786.__

Отпечатано в Отделе оперативной полиграфии химического факультета СПбГУ 198504, Санкт-Петербург, Старый Петергоф, Университетский пр., 26 Тел.: (812) 428-4043, 428-6919

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

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

04201358156

Райк Алексей Владимирович

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НАНОРАЗМЕРНЫХ СТРУКТУР НА ОСНОВЕ ОКСИДОВ МЕТАЛЛОВ

Специальность 05.13.18 — математическое моделирование, численные

методы и комплексы программ

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

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

Санкт-Петербург — 2013

Оглавление

Введение ....................................................................3

Глава 1. Математическая модель процессов адсорбции воды

на поверхности кристаллов........................................7

1.1 Уравнение Шредингера для многоэлектронных систем .... 7

1.2 Метод Хартри—Фока................................................9

1.3 Метод функционала электронной плотности ....................15

1.4 Выбор базисных функций..........................................25

1.5 Методы оптимизации геометрии кластерных систем............29

1.6 Кластерная модель расчета поверхностных свойств кристаллов MgO и ZnO......................................................34

Глава 2. Моделирование наноразмерных структур ..............36

2.1 Объекты исследования..............................................36

2.2 Кластерная модель поверхности (метод силового поля) .... 39

2.3 Расчеты на высокопроизводительном вычислительном комплексе в программном пакете Gaussian................................46

2.4 Кластерная модель поверхности (методы Хартри—Фока и функционала электронной плотности)..................................55

2.5 Взаимодействие кристаллических структур......................63

2.6 Модели полупроводников и диэлектриков на основе MgO

и ZnO................................................................66

Глава 3. Адсорбция на поверхности кристаллов MgO и ZnO . 70

3.1 Адсорбция молекулы воды на поверхности кристалла MgO . . 70 3.1.1 Адсорбция молекулы воды на поверхности кристалла

MgO (метод силового поля)................................73

СКу

3.1.2 Адсорбция молекулы воды на поверхности кристалла

(методы Хартри—Фока и функционала электронной плотности)..............................................75

3.2 Адсорбция молекулы воды на поверхности кристалла ZnO . . 86

3.3 Адсорбция молекулы воды на поверхности кристалла СаО . . 89

3.4 Адсорбция двух молекул воды на поверхности кристалла М^О 92 Глава 4. Расчет потенциалов взаимодействия молекулы воды

с поверхностью ......................................................97

4.1 Вычисление потенциалов взаимодействия........................97

4.2 Аппроксимация потенциалов аналитической функцией .... 103

4.3 Аппроксимация потенциалов с использованием численных методов .................................110

Заключение.................................116

Список литературы ...........................120

Приложение А. Сходимость процесса оптимизации геометрии 134

А1. Взаимодействие двух поверхностей................134

А2. Адсорбция молекулы воды над поверхностью..........134

АЗ. Две молекулы воды над поверхностью..............138

Приложение В. Директивы программного комплекса Gaussian 142

Введение

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

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

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

стижения поставленной цели были поставлены и решены следующие задачи:

1. Разработана методика моделирования кластерных фрагментов поверхности.

2. Рассчитаны энергетические характеристики поверхностных слоев.

3. Исследовано взаимодействие молекулы Н20 с поверхностью кристаллов MgO и ZnO.

4. Найдены потенциалы взаимодействия молекулы воды с поверхностями кристаллов.

Методами исследования являются методы математического и компьютерного моделирования, численного эксперимента, методы квантовой механики и физики твёрдого тела.

Положения, выносимые на защиту:

• Кластерные модели фрагментов поверхности с конечным числом атомов, исследуемые в программном пакете Gaussian.

• Способ расчета электронных характеристик, определяющих специфику поверхности кристаллов MgO и ZnO, а также их адсорбционные и каталитические свойства.

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

• Новая математическая модель механизма адсорбции молекулы воды на кристаллических поверхностях MgO и ZnO.

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

• Комплекс программ, реализующий созданную математическую модель.

Научная новизна. Впервые исследованы кластерные модели поверхностей, содержащие до 100 атомов поверхностного слоя. Найдены энергетические уровни кристаллических структур, определены значения зарядов на атомах, а также полная энергия систем.

Предложены новые структуры полупроводников с разными характеристиками запрещенной зоны на основе кубической кристаллической решетки 1^0 и ZnO с изменяемым процентным соотношением ионов магния и цинка.

Исследован механизм адсорбции молекул воды на различных кристаллических поверхностях ЇУ^О и ZnO. Показано, что при взаимодействии с поверхностью энергетически выгодным является образование водородной связи молекулы воды с атомом кислорода поверхности, тогда как адсорбция воды на поверхности кристалла 2пО обусловлена электростатическим взаимодействием.

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

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

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

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

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

Апробация работы. Результаты диссертационной работы докладывались и обсуждались на следующих конференциях: 40-я международная конференция студентов и аспирантов «Процессы управления и устойчивость» (СПб, СПбГУ, факультет ПМ-ПУ, 2009 г.); 41-я международная конференция студентов и аспирантов «Процессы управления и устойчивость» (СПб, СПбГУ, факультет ПМ-ПУ, 2010 г.); 43-я международная конференция студентов и аспирантов «Процессы управления и устойчивость» (СПб. СПбГУ, факультет ПМ-ПУ, 2012 г.); всероссийская конференция «Устойчивость н процессы управления» (КЗСР) им. В. И. Зубова (СПб, СПбГУ, 2010 г.); а также на научных семинарах кафедры моделирования электромеханических и компьютерных систем факультета прикладной математики - процессов управления СПбГУ.

Публикации. Основные результаты диссертационной работы отражены в 6 публикациях, список которых приведен в конце автореферата, в том числе 2 работы опубликованы в изданиях, рекомендованных ВАК Министерства образования и науки РФ.

Структура и объем работы. Диссертационная работа изложена на 133 страницах машинописного текста и состоит из введения, четырех глав, заключения и библиографического списка литературы, включающего 139 наименований. Работа содержит 31 рисунок, 19 таблиц и 2 приложения на 13 листах.

Глава 1. Математическая модель процессов адсорбции воды на поверхности кристаллов

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

1.1 Уравнение Шредингера для многоэлектронных систем

Основа теории электронной структуры вещества — это нерелятивистское стационарное уравнение Шредингера для многоэлектронной волновой

функции Ф:

#Ф(гь г2, • • •, Гдг, Къ К-2, • • •, Км) = ЯФ(гь г2, • ■ •, ГІУ, Иь И2, • • •, Н-м),

где Е — нолная энергия системы, Ф(гі, г2, ..., гдг, Ыь Яг, • • ■, К-м) — молекулярная волновая функция, (гь г2> • • ■) кь 1^2) • • • і Г^м) совокупность координат всех электронов и ядер, Й — оператор Гамильтона для молекулярной системы, состоящей из М ядер и N электронов при отсутствии магнитных и электрических полей. Гамильтониан Н имеет следующий вид:

n м

! 2 ^ МА

г=1 с*=1 л

n м „ n n м м „ „

ЕЕ^+ЕЕ^+ЕЕ^у- (ш)

і і ' га л . ' га „ 1ар

г=1 а=1 г=1 а= 1 /3>а

Здесь греческие индексы принадлежат ядрам, а латинские — электронам. Т и V являются операторами кинетической и потенциальной энергии соответственно. Оператор Т дополнен членом, учитывающим кинетическую энергию ядер, а оператор потенциальной энергии теперь содержит три слагаемых: энергию электрон-электронного отталкивания, энергию притяжения каждого электрона к каждому ядру молекулы и потенциальную энергию отталкивания положительно заряженных ядер.

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

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

Перепишем уравнение Шредпнгера в приближении Борна—Оппенгейме-

ра [1], в котором значительно более тяжелые ядра считаются неподвижными.

где г^- — положения электронов, 11а — положения ядер, ZCi — атомные номера ядер, /г, ше, е — фундаментальные константы, Е — полная энергия системы.

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

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

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

Е Ф = 0, (1.1.2)

1.2 Метод Хартри—Фока

онарнос состояние молекулы, которому соответствует энергия Ек- Электронная конфигурация молекулы, то есть совокупность ее молекулярных орби-талей, строится на основании двух фундаментальных принципов: принципа наименьшей энергии и принципа Паули. Это означает, что для описания электронной конфигурации молекулы с 2п или 2п — 1 электронами требуется п молекулярных орбиталей (или п энергетических уровней), так как на каждой орбитали могут располагаться два электрона [3]. Если орбитали вырождены (то есть соответствуют одному и тому же значению энергии), то они заполняются согласно первому правилу Хунда [4]. Орбиталь, па которой находятся два электрона, называется закрытой, в противном случае — открытой. Можно заметить, что концепция метода молекулярных орбиталей близка к концепции атомных орбиталей, что позволяет перенести математический аппарат, развитый для атомных орбиталей, па молекулярные орбитали. Отличие же этих методов заключается в том, что функции, отвечающие молекулярным орбиталям, являются многоцентровыми (по числу ядер в молекуле), а для атомных орбиталей они являются одноцентровыми.

Метод Хартри—Фока основан на использовании вариационного подхода [5]. Он позволяет, оставаясь в рамках приближения самосогласованного поля, поэтапно уточнять решения, расширяя используемый базис.

Запишем гамильтониан п-электронноп системы:

Английский математик Хартри в 1927 году предложил взаимодействие каждого электрона со всеми остальными заменить его взаимодействием с усредненным полем, создаваемым ядром и (п —1) электроном. Такой подход позволяет в (1.2.1) заменить потенциал е2/гкоторый зависит от координат двух электронов, выражением, описывающим межэлектронное взаимодействие как функцию координат каждого отдельного электрона. Согласно Хартри, многоэлектронную волновую функцию можно искать в виде простого произведения одноэлектронных волновых функций:

(1.2.1)

Ф(1, ...1п) = ф1(1)ф2{2)...'фп(п).

(1.2.2)

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