CC BY
95
Сер. 10. 2011. Вып. 2
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
ИНФОРМАТИКА
УДК 539.193
А. В. Райк, М. Е. Бедрина
МОДЕЛИРОВАНИЕ ПРОЦЕССА АДСОРБЦИИ ВОДЫ НА ПОВЕРХНОСТИ КРИСТАЛЛОВ
Введение. Вода практически всегда присутствует на поверхности твердых тел и влияет на ее свойства. Поэтому для практических целей важно теоретически исследовать взаимодействие воды с различными материалами. Процессы адсорбции воды на разных поверхностях изучаются давно и разнообразными методами [1-4]. Уникальные свойства воды, способность ее образовывать так называемые «водородные связи» приводят к сложной природе взаимодействия молекул воды с поверхностью. Поскольку образование водородной связи затрагивает перераспределение электронной плотности между частицами, корректное описание ее возможно только квантово-механическими методами. Важная информация о строении адсорбционных систем может быть получена на основе кластерных моделей [5-7], т. е. моделей, позволяющих представлять поверхность кристалла системой из конечного числа атомов.
Математическая модель. Для процесса взаимодействия адсорбированной молекулы с поверхностью математическая модель строится с привлечением аппарата квантовой механики. При этом удобно использовать кластерный подход при моделировании взаимодействия адсорбированных молекул с поверхностью благодаря простоте и гибкости модели, а также возможности применения стандартных методов квантовой химии.
Кластер, моделирующий фрагмент поверхности вместе с адсорбированной молекулой, рассматривается как единая замкнутая система, для которой учитываются все возможные виды взаимодействия: кинетическая энергия движения электронов и ядер как решение уравнения Лапласа, потенциальная энергия взаимодействия электронов с ядрами и друг с другом с учетом эффекта неразличимости электронов и их волновой природы. В адиабатическом приближении с разделением электронных и ядерных координат и представлением электронной волновой функции в виде линейной комбинации
Райк Алексей Владимирович — аспирант кафедры моделирования электромеханических и компьютерных систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Научный руководитель: доктор физико-математических наук,
проф. Н. В. Егоров. Количество опубликованных работ: 3. Научные направления: адсорбция на поверхности, математическое моделирование, численные методы. E-mail: raik.aleksey@gmail.com.
Бедрина Марина Евгеньевна — кандидат химических наук, доцент кафедры моделирования электромеханических и компьютерных систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 46. Научные направления: математическое моделирование, межмолекулярные взаимодействия. E-mail: m.bedrina@mail.ru.
© А.В. Райк, М.Е. Бедрина, 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 г.
