Моделирование образования химических связей при адсорбции

Логинова Мария Борисовна; Марамыгин Кирилл Вячеславович; Пономарев Александр Владимирович; Русинов Сергей Витальевич; Сакович Илона Олеговна

УДК 539.2

Моделирование образования химических связей при адсорбции

МГТУ им. Н.Э. Баумана, Калужский филиал, Калуга, 248000, Россия

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

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

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

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

1) первопринципные квантово-механические методы (ab initio);

2) классические методы молекулярной динамики (МД) на основе эмпирических потенциалов и законов классической физики;

3) полуклассические методы МД на основе квантово-механи-ческой аппроксимации сильных связей — методы МДСС.

Наиболее точными являются первопринципные квантово-меха-нические методы, базирующиеся на теории функционала плотности [1]. Однако вычислительные процедуры таких методов высоко-

затратны. Показателен рассмотренный в [2] пример поиска равновесной конфигурации аморфного кластера, содержащего 100 атомов. Для фиксированного набора внешних параметров при решении этой задачи методами МД потребовалось от 10 000 до 30 000 итерационных циклов, а время тестового расчета только одной итерации методами ab initio составило 8 ч.

Эффективные исследования большого круга задач межатомарных и молекулярных взаимодействий систем, содержащих более 100 элементов, как правило, проводят методами МД на основании хорошо апробированных эмпирических потенциалов. Стандартными потенциалами для анализа углеродных связей можно считать потенциалы Дж. Терсоффа [3, 4] и Д. Бреннера [5, 6]. Вместе с тем при использовании методов МД специфические особенности атомарных взаимодействий, обусловленные квантово-механической природой трансформации электронных структур, не поддаются воспроизведению и описанию, что в ряде случаев может приводить к неверным результатам. Разработанные в [7-9] методы МДСС существенно менее затратны по сравнению с методами ab initio и в то же время достаточно корректны, чтобы адекватно воспроизвести особенности изменения электронной плотности при образовании ковалентных связей. Детальное тестирование и анализ, проведенные в [8, 9], показали, что результаты, полученные на основе потенциалов сильной связи для различных структур и гибридизаций химических связей, находятся в очень хорошем соответствии как с известными экспериментальными данными, так и с результатами, полученными на основе методов ab initio.

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

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

rji = Rj + ri,

где Rj — координационный центр расчетной j-й ячейки; r — координата l-го атома в системе координат j-й ячейки.

Обозначим через hi (r) гамильтониан изолированного l-го атома, тогда волновое уравнение для атома с координатами rji примет вид

hl Ф ml (r - rjl) = Emi ф ml (r - rji ) ,

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

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

Н о = 2 (Г — ,),

, ,1

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

фтк =~1= 2 ехР (• к)фт1 (Г — Г,1), (1)

тШ

где N — число модельных ячеек; при этом собственные функции у ш гамильтониана Н записывают в виде линейной комбинации

У к = 2 Ст1Ф т1к . т1

Действуя на у к оператором Н (Н = Н0 + Нш), с учетом ортогональности функций Блоха получаем следующую систему линейных уравнений для вычисления собственных функций и собственных значений гамильтониана Н :

2 (Нт1 ,т7 ' — Ек^ тт'^ пп' ) Ст'I' (к ) = 0, (2)

т1

где Нт/,т1' представляет собой матрицу элементов (Фт/||Н|Фт7'к); Еш — собственные значения оператора Н. Подстановка волновых функций (1) в (2) позволяет получить:

N N

(-л " " ехр 1у,1 V гш / / \ I \\

к) = 22-N-\УтШ (Г — ,) Н Ут'1'к (Г — ,))

г (г, — г, п ' )•к

, г N

N

= 2 exp

j

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

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

H = Te + Üее + Üei,

где Te — кинетическая энергия электрона; Üee и Üei — потенциальная энергия взаимодействий электрон - электрон и электрон - ядро.

Таким образом, в соответствии с нотацией Хорсфилда уравнение Шрёдингера имеет вид

#| и) = е( n)\n),

где |n) и e(n) — собственные функции и собственные значения, при этом

и= 2 СИ и,

ia

здесь i, a — индексы соответственно сайта и орбиты.

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

Ей с (n) = с( n)c(n)

nia, jР^ jр _ ь ^ia 5

jP

в котором Hia jp — матричные элементы:

Ha, jp= (ia | H\ jp); (3)

2 CUxCk%m = 2 (n Iia) (ia Im) = 5 n,m.

ia ia

При ia ^ jp матричные элементы Hia,jp = (ia| H\ jp) представляют собой интегралы «перескока», в то время как Hia ia описывают

i (Rj + Г/ - Г/,)- k U m/k (Г - rj/ ) H V m'l k (Г - j')).

орбитальные энергии. В приближении сильной связи элементы матрицы (3) являются константами, которые находят с учетом следующих двух ограничений. Первое — энергетическое, оно связано с рассмотрением только тех атомарных орбит, для которых значения энергии электронов близки к энергии соответствующих валентных зон. (например, применительно к атомам углерода это орбиты 2s и 2p

(px, fy, pz)). Причем ненулевыми оказываются лишь такие интегралы «перескока», для которых угловые моменты относительно оси связи одинаковы. С учетом указанного обстоятельства только четыре интеграла «перескока» отличны от нуля: (ssa), (spa), (ppa) и (ррл:), где a соответствует орбитам с угловым моментом, равным нулю, а л — орбитам с угловым моментом ±1. Второе ограничение — дистанционное, оно связано с рассмотрением интегралов «перескока» только атомов, расстояние между которыми не превышает определенного радиуса обрезания.

Таким образом, с учетом рассмотренных ограничений недиагональные элементы матрицы Hia,jß = (ia| Й\ jß) (ia ^ jß) соответствуют устойчивым электронным топологиям для различных равновесных фаз кристаллических структур.

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

pia, jß=2 (4)

n

С учетом (4) выражение для энергии принимает вид

Etot = 2 2 P jß,iaHia, jß + Urep, ia,jß

тогда для сил, действующих на атомы, имеем

F dEtot J2V P DHia,jß ^DUrep

Fk =--:-=-i2 > p jß ia—г-+ ——

% %

Переходя непосредственно к моделированию перемещений взаимодействующих атомов средствами МД, целесообразно обратиться к следующей нотации. Результирующую компоненту силы Fia, действующую на i-й атом с координатами a (a = {x, y, z}), определим через градиент энергии и запишем классические уравнения движения в таком виде:

Л2т дЕ ^ = = ; (5)

Лг2 ду*

а

V- =-—

Уга л ,

здесь тга, уга — зависящие от времени г компоненты координат и скорости г-го атома массой т.

Систему уравнений (5) целесообразно решать на основе следующего корректирующего алгоритма. Полагая Г = тЯ, перепишем уравнение (5) в виде

Г = Я (г, г, г).

Первый шаг заключается в определении атомистических значений координат и скорости в момент времени I + Аг на основании атомистических значений координат и скорости в момент времени I -гМ, г = 0, к - 2, где к — порядок предикции.

Экстраполяцию атомистических координат и скоростей задают выражениями:

к-1

Г (г + Аг) = г (г) + Гг (г)Аг + Аг2 2 аЯ(г + [1 - г]М);

г=1 к-1

Г (г) Аг = г (г + Аг) - г (г) + Аг2 2 РгЯ(г + [1 - /]А/),

г=1

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

к-1 1

2 (1 - г) Д а г =-1-;

^ } г (д +1)(д + 2)' к 1 (6)

21 (1 - г)дРг = -^т, Д = 0,1-2.

г=1 Д + 2

Значения алгоритмических коэффициентов

Коэффициент к

1 2 3

а г 19 -10 3

Р г 27 -22 7

У г 3 10 -1

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

5г 7 6 -1

Коррекцию значений атомистических координат и скоростей проводят на основании выражений:

к-1

Г (Г + М) = г (Г) + Г (Г+ М2 $ У + [2 - 1]М);

г=1

к-1

(М )г (Г) = г « + М) - г () + А*2 $ ( + [2 - /]Д )(*),

г=1

где коэффициенты у , и 5, задают по аналогии с коэффициентами а, и Р^ — см. (6).

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

АЕ = Еа+в-( Еа + Ев). (7)

Необходимо обратить внимание на отрицательное, согласно определению (7), значение АЕ, причем чем более энергетически выгодна образующаяся связь, тем меньше, с учетом знака (АЕ < 0), абсолютное значение АЕ. Несмотря на энергетический выигрыш, образующийся в результате протекания химических реакций, в подавляющем большинстве случаев реагентам для инициации реакции требуется преодолеть определенный энергетический барьер Ебар, т. е. уже на этапе начального взаимодействия реагенты должны обладать достаточной избыточной энергией, превышающей значение Ебар.

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

ЕЛЯ = ЕЬО^ - .

Еа ~ Е% Е ;

= Ел? + Егз.

Е = еА + ЕВ - ЕШ ;

T->dis rebond ETS rebond I ETS . ETS \ . TTTS . Ea - Eg - E - Eg EA + EB j + E

■Jint ■>

eaas - efs - ets - ( ef+efbs )-( eas+ebs )+i

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

рисунок). Таким образом, первым (расчетным) этапом разрабатываемой методологии является определение уровней энергетических барьеров

Энергетическая диаграмма диссоциации адсорбируемого элемента на поверхности с учетом переходного (75) состояния (75, ¥5 — начальное и конечное состояния) [10]

( Ебар ).

возможных химических реакций в соответствии с иссле-

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

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

Второй этап связан непосредственно с моделированием кинематики взаимодействия элементов гетерогенной системы в рамках сформулированных на первом этапе необходимых термодинамических режимов и в интервалах установленных значений температуры (Т ...ти) и давления (Рн...Ри). Силы межатомарного взаимодействия расчитывают на

основе градиента полной энергии, потенциальную компоненту которой сначала определяют в МД-приближении, что позволяет решать кинематические уравнения движения для большого числа взаимодействующих элементов на основе алгоритма Верлета [11, 12]

с заданным шагом по времени (1.. .5 фс), а термодинамические условия, давление и температуру поддерживать на основе алгоритма Берендсена [12, 13]. Затем полученные в МД-приближении структуры уравновешивают на основе метода сильных связей [14], при этом для потенциальной энергии взаимодействия атомов находят матрицу плотности с помощью метода O(N) на основе оптимизации функционала третьего порядка [15].

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

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

ЛИТЕРАТУРА

[1] Paine M.C., Teter M.P., Allan D.C., Arias T.A., Joannopoulos J.D. Iterative Minimization Techniques for ab initio Total-Energy Calculations: Molecular Dynamics and Conjugate Gradients. Rev. Modern Phys, 1992, vol. 64, pp. 1045-1097.

[2] Patrice E.A., Gonis A., Colombo L. Tight-Binding Approach to Computational Materials Science. Boston, MIT, 1998, 542 p.

[3] Tersoff J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon. Phys. Rev. Lett, 1988, vol. 61, pp. 2879-2882.

[4] Tersoff J. Modeling Solid-State Chemistry: Interatomic Potentials for Multicomponent Systems. Phys. Rev. B, 1989, vol. 39, pp. 5565-5568.

[5] Brenner D.W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B, 1990, vol. 42, pp. 9458-9471.

[6] Brenner D.W., Shenderova O.A., Harrison J.A., Stuart S.J., Ni B., Sinnott S.B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys. Condens. Matter, 2002, vol. 14, pp. 783-802.

[7] Wang C.Z., Ho K.M. Tight-Binding Molecular Dynamics Studies of Covalent Systems. Adv. Chem. Phys, 1996, vol. XCIII, pp. 651-702.

[8] Wang C.Z., Ho K.M. Environment Dependent Tight-Binding Potential Model. Phys. Rev. B, 1996, vol. 53, pp. 979-987.

[9] Wang C.Z., Ho K.M. Material Simulations with Tight-Binding Molecular Dynamics. J. Phase Equil., 1997, vol. 18, pp. 516-527.

[10] Liu Z.P., Hu P. General Rules for Predicting where a Catalytic Reaction should Occur on Metal Surfaces. J. Am. Chem. Soc, 2003, vol. 125, pp. 1958-1967.

[11] Verlet L. Computer Experiments on Classical Fluids. Phys. Rev, 1968, vol. 159, pp. 98-103.

[12] Иванов В.А., Рабинович А.П., Хохлов А.Р. Методы компьютерного моделирования. Москва, Либроком, 2009, 662 с.

[13] Berendsen H.J.C., Postma J.P.M., Gunsteren W.F., Dinola A., Haak J.R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys., 1984. vol. 81, pp. 3684-3690.

[14] Xu C.H., Wang C.Z., Chan C.T., Ho K.M. A Transferable Tight-Binding Potential for Carbon. J. Phys. Condens. Matter, 1992, vol. 4, pp. 6047-6054.

[15] Li X.P., Nunes R.W., Vanderbilt D. Density-matrix Electronic-structure Method with Linear System-Size Scaling. Phys. Rev. B, 1993, vol. 47, pp. 10891-10894.

Статья поступила в редакцию 14.04.2015

Ссылку на эту статью просим оформлять следующим образом: Логинова М.Б., Марамыгин К.В., Пономарев А.В., Русинов С.В., Сакович И.О. Моделирование образования химических связей при адсорбции. Инженерный журнал: наука и инновации, 2015, вып. 4. URL: http://engjournal.ru/catalog/msm/nnm/1388.html

Логинова Мария Борисовна родилась в 1989 г., окончила Калужский филиал (КФ) МГТУ им. Н.Э. Баумана в 2011 г. Канд. физ.-мат. наук, доцент кафедры «Программное обеспечение ЭВМ, информационные технологии и прикладная математика» КФ МГТУ им. Н.Э. Баумана. Область научных интересов: физика конденсированного состояния, компьютерное моделирование и интеллектуальные системы.

Марамыгин Кирилл Вячеславович родился в 1986 г., окончил МГТУ им. Н.Э. Баумана в 2009 г. Канд. физ.-мат. наук, доцент кафедры «Программное обеспечение ЭВМ, информационные технологии и прикладная математика» Калужского филиала МГТУ им. Н.Э. Баумана. Область научных интересов: физика конденсированного состояния, компьютерное моделирование и проектирование программного обеспечения. e-mail: kirill.maramygin@gmail.com

Пономарев Александр Владимирович родился в 1989 г., окончил Калужский филиал (КФ) МГТУ им. Н.Э. Баумана в 2012 г. Аспирант кафедры «Программное обеспечение ЭВМ, информационные технологии и прикладная математика» КФ МГТУ им. Н.Э. Баумана. Область научных интересов: физика конденсированного состояния и компьютерные сети.

Русинов Сергей Витальевич родился в 1983 г., окончил Калужский филиал (КФ) МГТУ им. Н.Э. Баумана в 2006 г. Ассистент кафедры «Программное обеспечение ЭВМ, информационные технологии и прикладная математика» КФ МГТУ им. Н.Э. Баумана. Область научных интересов: физика конденсированного состояния, разработка и проектирование программного обеспечения. e-mail: rsage@mail.ru

Сакович Илона Олеговна родилась в 1992 г. Студентка кафедры «Программное обеспечение ЭВМ, информационные технологии, прикладная математика» Калужского филиала МГТУ им. Н.Э. Баумана. Область научных интересов: компьютерное моделирование, интеллектуальный анализ данных и распознавание образов. e-mail: ilona.sakovich@rambler.ru

Simulation of chemical bonds formation during adsorption

Bauman Moscow State Technical University, Kaluga Branch, Kaluga, 248000, Russia

The article describes main mathematical approaches to the description of interatomic interactions: quantum-mechanical methods, molecular dynamics methods based on empirical potentials and molecular dynamics methods based on approximation of tight-binding. Merits and demerits of the considered approaches are given. Theoretical part presents the mathematical apparatus for approximate calculation of the wave equation of the atom on the basis of the Hamiltonian and the theory of tight-binding. A correcting algorithm of calculation of kinematic characteristics is used for modeling movements of the interacting atoms. The authors consider two methodological stages: settlement on which initial conditions are set and levels of energy barriers of possible chemical reactions, and model stage where the assessment of forces of interatomic interaction is made are defined. A specific feature of modeling is in combination of the methods of molecular dynamics with approximation of tight-binding. The advantage of such approach is in significant increase of the size of the system without considerable loss of calculation accuracy. The methodology is demonstrated by an example of chemical bonds formation at adsorption.

Keywords: computer simulation, molecular dynamics methods, tight-binding approximation, heterogeneous systems

REFERENCES