Моделирование процесса адсорбции воды на поверхности кристаллов Текст научной статьи по специальности «Нанотехнологии»

Сер. 10. 2011. Вып. 2

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

ИНФОРМАТИКА

УДК 539.193

А. В. Райк, М. Е. Бедрина

МОДЕЛИРОВАНИЕ ПРОЦЕССА АДСОРБЦИИ ВОДЫ НА ПОВЕРХНОСТИ КРИСТАЛЛОВ

Введение. Вода практически всегда присутствует на поверхности твердых тел и влияет на ее свойства. Поэтому для практических целей важно теоретически исследовать взаимодействие воды с различными материалами. Процессы адсорбции воды на разных поверхностях изучаются давно и разнообразными методами [1-4]. Уникальные свойства воды, способность ее образовывать так называемые «водородные связи» приводят к сложной природе взаимодействия молекул воды с поверхностью. Поскольку образование водородной связи затрагивает перераспределение электронной плотности между частицами, корректное описание ее возможно только квантово-механическими методами. Важная информация о строении адсорбционных систем может быть получена на основе кластерных моделей [5-7], т. е. моделей, позволяющих представлять поверхность кристалла системой из конечного числа атомов.

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

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

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

проф. Н. В. Егоров. Количество опубликованных работ: 3. Научные направления: адсорбция на поверхности, математическое моделирование, численные методы. E-mail: [email protected].

Бедрина Марина Евгеньевна — кандидат химических наук, доцент кафедры моделирования электромеханических и компьютерных систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 46. Научные направления: математическое моделирование, межмолекулярные взаимодействия. E-mail: [email protected].

© А.В. Райк, М.Е. Бедрина, 2011

по заданным базисным функциям (атомным орбиталям) решается уравнение Шредин-гера

НФ (г, т) = ЕФ (г, т),

где Е - полная энергия системы; Ф(г,т) - молекулярная волновая функция; Н - оператор Гамильтона для молекулярной системы.

Гамильтониан Н имеет следующий вид:

1 1 М 1

й = t + v = — Vv2--V —v д -

^ 2 2 ^ МА А

г=1 а=1

N М N N М М

-EE^+EEè+EEff-

- 1 1 'га . ... ' гп л Rati

г=1 а=1 г=1 п>г J a=1 ¡3>а

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

Опираясь на основную идею метода функционала электронной плотности (DFT) [8], что все свойства атомов и молекул определяются распределением электронной плотности в системах, запишем отдельные компоненты полной энергии основного состояния:

Е0(р0) = Т (ро) + Еее(ро) + ENe(p0)ï (1)

здесь

p(r) = N J■■■ J |^(xb x2, ..., xn)|2 ds1 dx2 ... dxN

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

Разделим выражение (1) на две части. В первую входит потенциальная энергия электрон-ядерного притяжения ENe(po) = J Po(r)VNe dr, как часть, зависящая от конкретной системы; во вторую — Т(ро) и Eee(po) как универсальные, т. е. не зависящие от N (число электронов), Ra (позиции ядер) и Za (заряды ядер). В результате получаем

Е0(Р0) = j P0(r)VNe dr + Т(Р0) + Eee(P0).

Объединив независимые части в новую переменную Fhk(P0), приходим к

Е0Р) = J P0(r)VNe dr + Fhk(P0).

Функционал Eee(р) можно представить в виде суммы двух функционалов:

Еее(р) = Ô /У Р<У ^Р<У ^ dr 1 dr2 + ЕпсЛ(р) = J(p) + ЕпгЛ(р),

¿JJ '12

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

Выпишем функционал Енк

Енк(р) = Тз (р) + 3 (р)+Ехс (р). (2)

В (2)

1 "

тв = ~^№ 1у21^>>

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

Ехс(р) = (Т(р) - Тз(р)) + (Еее(р) - 3(р)) = Тс(р) + Епс1(р).

Дополнительная часть функционала кинетической энергии Тс представляет собой поправочный член.

Тогда энергию основного состояния можно записать в виде

Е (р(г)) = Тз (р) + J (р) + Ехс (р) + ENe(p) =

1 М 1 М М ГГ 1

= II №(Г1)|2—№(Г2)|2*1*2 +

і 3

А

Далее, применив вариационный принцип, получим результирующие уравнения

4У2 +

Р(Г2) л , 17 / \ А

------(ІГ2 + ^хс(гі) - > ^---------

Г12 ^г' Г1А

фт —

= | ^ ^ ^ і Фт = ЄтФті

где

^2^т = Те (р) + I Ує«(г)р(г) ¿г.

Выбор формы Ехс напрямую влияет на точность приближения. В настоящем исследовании использовали гибридный функционал, предложенный в работах [9, 10]. Волновая функция Ф(г, т) представлялась суперпозицией гауссовых функций, 6 из которых - атомные орбитали остовных электронов, а 3 и 1 - атомные орбитали валентных электронов.

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

Методика построения кластерной модели. В процессе отработки построения кластерной поверхности были проведены расчеты моделей с разным количеством атомов в поверхностном слое и различным количеством слоев. В частности, рассматривались кластерные структуры размером [3x3], [5x5], [7x7] (одно- и двухслойные) и [9x9] (одно-, двух-, трех- и четырехслойные). Расчеты проводились методами ИНТ/6-3Ю, B3LYP/6-31G с применением программного пакета СаиБ81ап-03 [11] на высокопроизводительном вычислительном комплексе факультета прикладной математики-процессов управления СПбГУ.

В качестве начальной геометрии выбиралась симметричная и слегка искаженная (полученная методом молекулярной механики ММ+) структура поверхности с различными межатомными расстояниями 2.0 и 1.77 А. Для всех моделей была определена равновесная геометрия и найдена полная энергия системы [12]. Независимо от выбранной начальной геометрии, результатом оптимизации являлась поверхность со средними межатомными расстояниями порядка 2.0—2.1 А.

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

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

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

Таблица 1. Энергии кластерных моделей Е и Е/п

Слой Е, а. е. Е/п, а. е.

1 • [3 X 3] 1 301.5539 144.6171

2 • [3 X 3] 2 478.5192 137.6955

1 • [5 х 5] 3 504.9294 140.1972

2 • [5 X 5] 6 885.8826 137.7177

1 • [7 X 7] 6 685.1640 136.4319

2 • [7 X 7] 13 496.9676 137.7241

Адсорбция молекулы воды на поверхности кристаллов. Модель с молекулой Н2О над поверхностью кристалла представляет собой супермолекулярную систему, в которой учитываются взаимодействия всех атомов и конкурируют процессы отталкивания, притяжения, короткодействующие эффекты обмена электронов. В этом случае матрица координат задавалась двумя способами: в виде Z-матрицы (матрицы внутренних координат), позволяющей фиксировать координаты атомов, входящих в состав кластера, и в декартовых координатах. При рассмотрении взаимодействия молекулы воды с поверхностью, которая задается «фиксированной» геометрией, молекула воды соединяется с атомами кластера с помощью соединительного атома, позволяющего сохранить целостность структуры Z-матрицы, никак не влияющего на нахождение

равновесной геометрии. Полная оптимизация структуры кластера в процессе адсорбции Н2О приводит к искажению решетки и потере симметрии.

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

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

ДЕсвязи — Еполн (Епов + ЕЯ2 О \

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

Модель 2-[7x7] поверхности кристалла MgO с адсорбированной молекулой воды

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

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

Водородная связь - это разновидность донорно-акцепторной связи, т. е. невалентное взаимодействие, которое возникает между двумя атомами кислорода через водород [13].

Значения вычисленных полных энергий системы (М^0)крист-Н20 также во всех случаях свидетельствует о более выгодном энергетическом состоянии образования водородной связи молекулой воды с атомом кислорода поверхности. По-видимому, выявленное поведение молекулы воды над поверхностью можно объяснить экранированием ионов магния делокализованными р-орбиталями кислорода. Кажущееся очевидным электростатическое взаимодействие дипольной молекулы воды с ионами магния оказывается менее выгодным. Следовательно, корректное рассмотрение короткодействующих взаимодействий в модели адсорбции воды над поверхностью кристалла необходимо проводить с учетом обменных эффектов и электронной корреляции, которая косвенно учитывается в расчетах при использовании гибридного потенциала B3LYP.

Как показали проведенные расчеты, выбор начальной точки расположения молекулы воды над поверхностью кристалла не влияет на результат определяемой равновесной геометрии. Молекула воды в равновесном, наиболее энергетически выгодном положении образует водородную связь с атомом кислорода поверхности М^О, при этом молекула воды в поле поверхности незначительно искажается. Частоты колебаний в этом случае заметно изменяются.

Данные об изменении молекулярной структуры воды при образовании водородной связи приведены в табл. 2.

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

Параметр О изолированная 1 12 ^ и;| Д поверхностью

Д(М^..О), А 2.3705

д(оп ...но, А 1.69555

д(нь. .о), А 0.9759 1.02015

д(о-н2), А 0.9759 0.97383

¿(Щ-О-Нз) 108.304 111.502

9н, , а. е. 0.359 0.525

до, ^ е. -0.718 -0.750

<?Н2, а. е. 0.359 0.415

Анализ потенциальной поверхности системы (MgO)крист-H2O показывает, что атом водорода колеблется по линии водородной связи между двумя энергетическими минимумами, при этом расстояние в молекуле воды Н]—Оп изменяется от 1.1 до 1.5 А, что предполагает возможность диссоциации адсорбированной молекулы воды через разрыв ОНх-связи.

При изучении адсорбции молекулы воды на поверхности кристалла ZnO (100) применялась аналогичная методика расчетов. В качестве модели поверхности выбиралась двухслойная кластерная модель [7x7]. Расстояния Zn-O, полученные в результате предварительной оптимизации геометрии кластера (1-[7х7]), составляют 1.96-1.97 А. Процесс адсорбции рассматривался над «замороженной» поверхностью с вышеуказанными геометрическими значениями. Анализ движения молекулы по потенциальной кривой свидетельствует о существенном отличии от случая

с поверхностью М^О. Для кристалла оксида цинка не наблюдалось стремления к образованию водородной связи с кислородом поверхности, как в кристалле MgO. То есть поверхность кристалла является более гидрофобной, это косвенно подтверждается и экспериментальными данными. Природа взаимодействия в данном случае носит электростатический характер. Если сравнить распределение величин зарядов, полученных на поверхности для кристаллов MgO (1.728 а.е. для атома магния, —1.163 а.е. для атома кислорода) и ZnO (1.432 а.е для атома цинка, —0.956 а.е. для атома кислорода), можно видеть, что величины зарядов Mg и О в кристалле MgO больше, чем в ZnO (можно говорить, что в MgO связь более ионная). В целом в случае кристалла ZnO молекула воды располагается в соответствии с электростатическим принципом: атомом кислорода над атомом цинка, а атомы водорода ориентированы к атомам кислорода поверхности (но находятся на значительном расстоянии). В табл. 3 приведены данные об изменении геометрических параметров адсорбированной молекулы воды.

Таблица 3. Рассчитанные межъядерные расстояния и заряды на атомах для адсорбированной молекулы воды над поверхностью кристалла ZnO

Параметр Н2 Оизолированная Н20над поверхностью

Я(Ні-О), А 0.9759 0.9807

Д(о-н2), А 0.9759 0.9706

¿(Ні-О-Нз) 108.304 108.254

дщ, а. е. 0.359 0.445

ЯО, а. е. -0.718 -0.680

9Н2 > а. е. 0.359 0.439

Сравнение с экспериментом. Критерием достоверности предлагаемой модели является сравнение вычисляемых с ее помощью характеристик системы с экспериментальными данными и с данными, полученными в других модельных подходах. Прочность связи молекулы воды с поверхностью для различных «замороженных» молекулярных кластеров иллюстрирует табл. 4. Приведенные значения хорошо совпадают с эмпирической оценкой энергии связи, выполненной методом корреляций (14— 17 ккал/моль), так как прямые экспериментальные измерения энергии связи крайне сложны [14].

Таблица 4. Энергии связи молекулы H2O с поверхностью кластеров MgO

Слой ®п+Н20, ккал/моль ДЕ, ккал/моль

1 • [3 X 3] 1 377.9638 14.9709

2 • [3 X 3] 2 554.9277 14.0617

1 • [5 X 5] 3 504.9294 18.6683

2 • [5 X 5] 6 962.3020 20.9415

1 • [7 X 7] 6 886.6838 22.3315

2 • [7 X 7] 13 571.3390 26.5496

В табл. 5 вычисленные значения частот нормальных колебаний свободной и адсорбированной молекулы воды сравниваются с данными ИК-спектроскопии [15], СХПЭЭ-ВР (спектроскопия высокого разрешения характеристических потерь энергии электронами) [16] и расчетами кристаллической структуры в базисе плоских волн [17]. Возможно, экспериментально определяемая частота 3369 см-1 является суперпозицией симметричного ^2 и антисимметричного VI валентных колебаний ОН-связей, тем

самым среднее значение вычисленных частот будет находиться еще ближе к экспериментальному. Полученная авторами работы [17] частота 3850 см-1 существенно завышена по сравнению с экспериментальной.

Таблица 5. Частоты нормальных колебаний (см 1)

Параметр H2 Оизолированная Н2С)Над поверхностью

v\ VI эксперимент 3781 3756 [15] 3738

V2 ^эксперимент 3616 3657 [15] 3138 3369 [16] 3850 [17]

Vi ^эксперимент 1619 1595 [15] 1611 1650 [16]

Заключение. Для понимания процессов, происходящих на поверхности, необходима информация о механизме взаимодействия адсорбированных молекул с кристаллом. В данной работе разработана методика моделирования кластеров и изучено взаимодействие молекул воды с поверхностью. Было показано, что в качестве расчетной модели достаточно использовать двухслойный кластер [7x7]. Приведено сравнительное рассмотрение взаимодействия молекул воды с поверхностями кристаллов MgO и ZnO. Результаты исследований свидетельствуют о различной природе взаимодействия. С поверхностью (100) MgO молекула воды образует устойчивую водородную связь. В случае поверхности ZnO этого не происходит.

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

Литература

1. Абаренков И. В., Третьяк В. М., Тулуб А. В. О механизме диссоциации молекул воды на (100) поверхности кристалла MgO в модели парных неэмпирических потенциалов // Хим. физика. 1985. Т. 4, № 4. С. 974-980.

2. Bandura A. V., Sykes D. G. Adsorption of Water on the Ti02 (Rutile) (110) Surface: A Comparison of Periodic and Embedded Cluster Calculations // J. Phys. Chem. B. 2004. Vol. 108. P. 7844-7853.

3. Эварестов Р. А., Бандура А. В. Компьютерное моделирование адсорбции воды на поверхности кристаллических оксидов титана, олова, циркония и гафния // Рос. хим. журн. 2007. Т. 5. C. 149-158.

4. Fox H., Gillan M. J., Horsfield A. P. Methods for calculating the desorption rate of an isolated molecule from a surface: Water on Mg0(001) // Surface Science. 2007. Vol. 601. P. 5016-5025.

5. Bellert D, Burns K. L., Wampler R. et al. An accurate determination of the ionization energy of the MgO molecule // Chem. Phys. Lett. 2000. Vol. 322. P. 41-44.

6. Puente de la E., Aguado A., Ayucla A. et al. Structural and electronic properties of small neutral

(Mg0)n clusters // Phys. Rev. B. 1997. Vol. 56. P. 7607-7614.

7. Malliavin M. J., Coudray C. Ab initio calculations on (Mg0)n, (Ca0)n, and (NaCl)n clusters (n =

1-6) // J. Chem. Phys. 1997. Vol. 106. P. 2323-2329.

8. Koch W., Holthausen M. C. A Chemist’s Guide to Density Functional Theory. Ed. 2. Weinheim: Wiley-VCH, 2002. 293 p.

9. Becke A. D. Density-functional thermochemistry. III. The role of exact exchange // J. Chem. Phys. 1993. Vol. 98. P. 5648-5652.

10. Lee C., Yang W., Parr R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density // Phys. Rev. B37. 1988. P. 785-789.

11. Frisch M. J., Trucks G. W., Schlegel H. B. et al. GAUSSIAN-03, Rev. B.05. Pittsburgh PA: Gaussian, 2003. 596 p.

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

13. Каплан И. Г. Введение в теорию межмолекулярных взаимодействий. М.: Наука, 1982. 312 с.

14. Henrich V. E., Cox P. A. The Surface Science of Metal Oxides. Cambridge, England: Cambridge University Press, 1994. 461 p.

15. Benedict W. S., Gailar N., Plyler E. K. Rotation-Vibration Spectra of Deuterated Water Vapor // J. Chem. Phys. 1956. Vol. 24. P. 1139-1165.

16. Wu 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.

17. Wang Y., Nguyen H. N., Truong T. N. Mechanisms of and Effect of Coadsorption on Water Dissociation on an Oxygen Vacancy of the MgO(100) Surface // Chem. Eur. J. 2006. Vol. 12. P. 58595867.

Статья рекомендована к печати проф. Н. В. Егоровым.

Статья принята к печати 16 декабря 2010 г.