Научная статья на тему 'Моделирование потенциалов межмолекулярного взаимодействия'

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

CC BY
395
77
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
КВАНТОВАЯ МЕХАНИКА / ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ / ПОТЕНЦИАЛЫ МЕЖМОЛЕКУЛЯРНОГО ВЗАИМОДЕЙСТВИЯ / АДСОРБЦИЯ / DFT / QUANTUM MECHANICS / PARALLEL COMPUTING / INTERMOLECULAR INTERACTION POTENTIALS / ADSORPTION

Аннотация научной статьи по нанотехнологиям, автор научной работы — Райк Алексей Владимирович, Егоров Николай Васильевич, Бедрина Марина Евгеньевна

Вычислены потенциалы межмолекулярного взаимодействия молекулы воды с поверхностью (100) кристаллов MgO и ZnO. В качестве расчетных методов исследования применялись квантово-механический метод B3LYP/6-31G, косвенно учитывающий корреляцию электронов, метод атом-атомных потенциалов. Потенциалы взаимодействия воды с поверхностью аппроксимированы аналитической функцией Букингема. Проведено сравнение результатов, полученных с помощью аналитических расчетов и применения нелинейного метода наименьших квадратов Левенберга-Марквардта. Аналитический вид потенциалов может быть использован в статистических методах Монте-Карло и молекулярной динамики для исследования распределения большого числа молекул над поверхностью с учетом изменения температуры.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по нанотехнологиям , автор научной работы — Райк Алексей Владимирович, Егоров Николай Васильевич, Бедрина Марина Евгеньевна

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

Modelling intermolecular interaction potentials

Potentials of intermolecular interaction of a molecule of water with MgO and ZnO crystals (100) surface were calculated. The B3LYP/6-31G quantum mechanics method, which takes into account electron correlation, and the atom-atom potential method were used for calculations. Potentials of the water molecule interaction with the surface were approximated by the Buckingham analytic function. Analytical calculations were compared with the results obtained using the LevenbergMarquardt method for non-linear least squares. Analytical expressions of potentials can be used to carry out statistical calculations using Monte-Carlo and molecular dynamics methods, which allow analyzing large numbers of molecules at a surface taking into account temperature changes.

Текст научной работы на тему «Моделирование потенциалов межмолекулярного взаимодействия»

Сер. 10. 2012. Вып. 3

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

ИНФОРМАТИКА

УДК 539.193

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

МОДЕЛИРОВАНИЕ ПОТЕНЦИАЛОВ МЕЖМОЛЕКУЛЯРНОГО ВЗАИМОДЕЙСТВИЯ

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

В настоящей работе были рассмотрены потенциалы межмолекулярного взаимодействия молекулы воды с поверхностью (100) кристаллов MgO и ZnO. В качестве методов исследования применялись методы атом-атомных потенциалов и квантово-механический с учетом корреляции электронов [2].

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

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

Для наиболее энергетически выгодной и актуальной в практическом смысле ориентации (а) были также вычислены квантово-механические потенциалы взаимодействия H2O с поверхностями MgO и ZnO. Расчеты этих потенциалов проводились методом функционала электронной плотности с гибридным потенциалом B3LYP/6-31G

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

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

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

© А. В. Райк, Н.В. Егоров, М.Е. Бедрина, 2012

Е, а.е.

9436

2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 5.8 6

г, А

Рис. 1. Графики потенциалов взаимодействия, полученные для различных ориентаций

молекулы воды над поверхностью кристалла MgO Ориентация молекулы воды: 1 - оси горизонтально (случай а), 2 - оси вертикально (случай Ь), 3 - в вертикальном положении (случай с).

в программе СаиББ1ап-03 [3] на высокопроизводительном вычислительном комплексе факультета прикладной математики-процессов управления СПбГУ [4, 5]. Для установления значения полной энергии Е в точке Я^ требуется около 19 ч машинного времени. Следует отметить, что при этом подходе не рассматривалось образование водородной связи молекулой воды с поверхностью кристалла MgO, что, как было показано ранее [5], является наиболее вероятным.

При определении положения адсорбированной молекулы над двухслойной поверхностью оксидов магния и цинка молекула воды в начальной точке располагалась на расстоянии 2 А над поверхностью. В равновесной конфигурации расстояние г Mgп-Он2о равно 2.9 А, а г Znп-Он2о составляет 3.4 А, что соответствует точкам минимума энергии рассчитанных потенциалов.

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

= Се— - § - ^

Я*

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

Потенциал Букингема основан на предположении об экспоненциальной зависимости сил отталкивания между молекулами от расстояния между ними. Он представляет собой сумму трех слагаемых. Первое слагаемое, Св-пп, характеризует отталкивание, диполь-дипольное взаимодеиствие, третье, - диполь-квадрупольное

второе,

в

д6'

взаимодействие [6]. По сравнению с потенциалом Леннарда-Джонса эта форма более сложна для математической обработки в связи с наличием одновременно экспоненциальной и степенной зависимостей, в то же время она более реалистична физически.

Для того чтобы аппроксимировать найденные потенциалы, нужно определить неизвестные константы при известных значениях К и у>(К). Так как имеется четыре неизвестные константы, то была составлена система из четырех нелинейных уравнений с четырьмя неизвестными и решена одним из доступных способов:

р. А

= дё-дё, (1)

Р А

(р(Ц2) = Се-1}1Ъ- — (2)

Щ Щ'

рА

^(Д 3) = Се-™з____ (3)

рА

<рт = Се-°в+- — -— (4)

Из уравнения (1)

Кл Кл

рА

Кб

константа А выражалась через неизвестные константы Р, С, Б и известные параметры К? и ц>(К?):

Полученное выражение для А подставлялось в уравнение (2):

-(^(ДО + ^-Се-^Д?

Р

-156 =ср(Н2),

г?8 рб

К2 К2

в котором константа Р выражалась через неизвестные С, Б и известные К?, К2, ^(К?), Д| [Се^ - ^(Д2)] - Д? [Се-73«! - *>№)]

о =

К2 - К2

В уравнении (3)

К|К2 - К?К2

(к? (К2 - К2) [Се-™1 - ^(К1)] +

+ К? (К2 - К?) [Се-^ - + СК? (К2 - К22) = ^(К3)

константа С выражалась через неизвестный параметр Б и известные константы:

1

С

К? (К2 - К2) е-ОК1 + К? (К2 - К2) е-сд2 + К? (К? - К2) е-СДз

х К?^(К!) (К2 - К2) + К?^(К2) (К2 - К2) + К|^(К3) (К2 - к2)

и подставлялась в (4):

1

я? (Я? - Я? ) в-оя1 + Я? (Я2 - Я? ) в-оя2 + я?? Я? (Я? - Я? ) е-^з

х (-Я? [Я?^(Я?) (Я? - Я?) + Я?V (Яз) (Я? - Я?) "Я?V (Я!) (Я? - Я?) + Я?V (Яз) (Я? - Я?)" Я? V (Я!) (Я? - Я?) + Я?V (Я2) (Я? - Я?)

+ Я? - Я?

е-ЯД1 +

а-ОН2

-ОЯз

+

+ Я?

Я?V (Я!) (Я? - Я?) + Я?V (Я?) (Я? - Я?) +

+ Я?V (Яз) (Я? - Я?)

-оял _

= V(Я4 ),

которое представляет собой уравнение с одним неизвестным параметром Б. Решая это уравнение, можно получить значение для Б и подставить его в (3). Теперь уравнение (3) является уравнением с одним неизвестным С. Константа В находится из уравнения (2) после подстановки в него уже определенных С и Б. И, наконец, подставляя В, С и Б в (1), можно получить значение для А, а вместе с ним и вид аппроксимирующей функции.

9436

2.3 2.5 2.7 2.9 3.1 3.3 3.5 3.7 3.9 4.1 4.3 4.5 4.7 4.9 5.1 5.3 5.5 5.7 5.9 6.1 6.3

г, А

Рис. 2. Аппроксимация потенциала MgO аналитической функцией (1) 2 - численный эксперимент.

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

СВ

А

в

Здесь дополнительный член -^то характеризует квадруполь-квадруиольное взаимодействие (вносит небольшой вклад в изменение потенциала). Дополнительное слагаемое предоставляет еще одну точку на графике, через которую будет проходить аппроксимирующая функция. Наличие дополнительной точки, с одной стороны, является преимуществом, так как позволяет более точно приблизиться к аппроксимируемой функции, с другой - значительно усложняет задачу поиска вида аппроксимирующей функции. В этом случае для нахождения пяти неизвестных констант решается система из пяти нелинейных уравнений. Такая система была решена способом, аналогичным вышеизложенному, и были получены аналитические значения для констант А, В, С, Б и Е.

В рамках этой задачи было проведено сравнение квантово-механических потенциалов для систем MgO-H2О и ZnO-H2О (рис. 3).

Рис. 3. Сравнение потенциалов, полученных методом B3LYP/6-31G для систем И2С^пС (1) и ИС^С (2)

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

Аппроксимация потенциалов с использованием численных методов. С помощью программных пакетов Gnuplot и Origin была произведена численная аппроксимация полученных значений энергии функцией Букингема. В этих пакетах реализован нелинейный метод наименьших квадратов Маркварта-Левенберга, представляющий собой комбинацию метода наискорейшего спуска по градиенту и метода Гаусса-Ньютона. Пусть имеется задача о наименьших квадратах вида

m m

F (x) = fi (x) = (x) - f (x))2 ^ min, x = (xi x, ...,xn).

i=1 i=1

Выпишем градиент

УГ(х) _ 27т(х)/(х).

д/(хз)

Здесь </(х) - матрица Якоби вектор-функции Г(х), т. е. ./(х) = —--, где 1 ^ ] ^ то,

ОХ1

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

1 ^ г ^ п. Матрица Гессе для компоненты /(х) будет иметь следующий вид:

Н(х) _ 7т(х)7(х).

Тогда, согласно методу Гаусса-Ньютона, очередное направление р находится из системы

7т(х)7(х)р _ -7т(х)/(х).

Направление поиска Марквардта-Левенберга определяется из системы

(Зт(хк) 7(хк) + ХкI) ри _ -7Т(хк)/(хк),

где Хк - некоторая неотрицательная константа, своя для каждого шага; I - единичная матрица.

Выбор Хк можно осуществлять, делая его достаточным для монотонного спуска Г(х), т. е. увеличивать параметр до тех пор, пока не будет достигнуто условие Г (хи+!) < Г (хк).

Нетрудно заметить, что при Хи _ 0 алгоритм вырождается в метод Гаусса-Ньютона, а при достаточно большом Хк направление рк незначительно отличается от направления наискорейшего спуска. Следовательно, при правильном подборе параметра Хи можно добиться монотонного убывания минимизируемой функции. Неравенство Г(хи+1) < Г(хи) всегда можно обеспечить, выбрав Хи достаточно большим. При этом теряется информация о кривизне, заключенная в первом слагаемом, и проявляются все недостатки метода градиентного спуска: в местах пологого наклона антиградиент мал, а в местах с крутым наклоном - велик, в то время как в первом случае желательно делать большие шаги, а во втором - малые. Для учета информации о кривизне единичную матрицу можно заменить на диагональ матрицы Гессе. Таким образом можно достичь увеличения шага вдоль пологих участков и уменьшения вдоль крутых спусков:

(7т(хи(хи) + Хиdiag 7т(хи(хи)]) ри _ -7т(хи)/(хи). (5)

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

Таблица 1. Значения констант А, В, С, Б для двух систем

Константа Gnuplot Origin

Система, Hg О ■ МдО

А -1.75559 ± 0.2358 -1.75481 ± 0.23618

В -1.25938 ± 0.1237 -1.25981 ± 0.124

С -0.209011 ± 0.01215 -0.20906 ± 0.01223

D 0.812155 ± 0.01694 0.81223 ± 0.01714

Система, О ■ ZnO

А -1.39858 ± 0.1739 -1.39271 ± 0.15080

В -2.13085 ± 0.06912 -2.13085 ± 0.06841

С -0.230875 ± 0.004071 -0.22721 ± 0.00911

D 0.806972 ± 0.00471 0.79913 ± 0.00625

Из табл. 1, 2 видно, что значения А, В, С и Б в целом совпадают. Но количество

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

Координата Потенциал MgO Остатки MgO Потенциал ZnO Остатки ZnO

1.5 0.1163227 -0.0251267

1.6 0.0605802 -0.0163898

1.7 0.0257222 -0.0125163

1.8 0.0043342 -0.0102616 0.0213224 -0.0000878

1.9 -0.0086180 -0.0086803 0.0036973 0.0001588

2 -0.0159054 -0.0070708 -0.0072097 0.0000844

2.1 -0.0197246 -0.0055513 -0.0138611 -0.0000632

2.2 -0.0213742 -0.0041429 -0.0177737 -0.0000952

2.3 -0.0217497 -0.0029377 -0.0199036 -0.0001102

2.4 -0.0213818 -0.0019510 -0.0208654 -0.0000683

2.5 -0.0206117 -0.0011880 -0.0210611 -0.0000150

2.6 -0.0196346 -0.0006213 -0.0207576 0.0000394

2.7 -0.0185809 -0.0002323 -0.0201339 0.0000686

2.8 -0.0175244 0.0000056 -0.0193114 0.0000759

2.9 -0.0165056 0.0001200 -0.0183729 0.0000667

3 -0.0154743 0.0002076 -0.0173748 0.0000454

3.1 -0.0146199 0.0001103 -0.0163559 0.0000269

3.2 -0.0137949 -0.0000033 -0.0153422 0.0000256

3.3 -0.0129829 -0.0001030 -0.0143514 0.0000270

3.4 -0.0122029 -0.0001987 -0.0133949 -0.0000366

3.5 -0.0114722 -0.0003026 -0.0124798 -0.0000305

3.6 -0.0107167 -0.0003374 -0.0116104 -0.0000247

3.7 -0.0100569 -0.0004226 -0.0107887 -0.0000247

3.8 -0.0094363 -0.0005016 -0.0100153 -0.0000279

3.9 -0.0088412 -0.0005616 -0.0092896 -0.0000318

4 -0.0082655 -0.0005979 -0.0086106 -0.0000330

4.1 -0.0077157 -0.0006187 -0.0079765 -0.0000284

4.2 -0.0071952 -0.0006293 -0.0073854 -0.0000171

4.3 -0.0066998 -0.0006277 -0.0068352 0.0000030

4.4 -0.0062289 -0.0006154 -0.0063237 0.0000355

4.5 -0.0057769 -0.0005888

4.6 -0.0053419 -0.0005483

4.7 -0.0049231 -0.0004949

4.8 -0.0045217 -0.0004319

4.9 -0.0041375 -0.0003609

5 -0.0037700 -0.0002831

5.1 -0.0034148 -0.0001959

5.2 -0.0030780 -0.0001068

5.3 -0.0027595 -0.0000173

5.4 -0.0024612 0.0000695

5.5 -0.0021831 0.0001521

5.6 -0.0019248 0.0002299

5.7 -0.0016855 0.0003025

5.8 -0.0014650 0.0003691

5.9 -0.0012626 0.0004293

6 -0.0010786 0.0004821

итераций, которое потребовалось программным пакетам для получения подобных результатов, существенно отличается. Для пакета Gnuplot оно колеблется в пределах 20-30, для пакета Origin равно 200-250. Такие отличия можно объяснить различными компьютерными реализациями алгоритма Марквардта-Левенберга. Следует отметить, что данные наборы констант являются не единственными. Итерационный процесс сходится к одному из минимумов. На конечный результат большое влияние оказывают

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

А _ В _ С _ -1, Б _1. При начальных значениях А _ В _ С _ Б _10 итерационный процесс сходится к другому минимуму, однако величины констант не соответствуют физическому смыслу слагаемых потенциала Букингема.

а.е. а

............................................

2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 ■ 1 1 4 4.2 4.4 4.6 4.8 5 5.2

. 1 1

" 1 -

\

- 1 ¿А

' \ У

\ Я К Л-

^ лГ

1.8 \2 2.2 2.4 2.6 2.8

3.2 3.4 3.6 3.8

4 4.2 4.4 4.6

г Л

...о.. 7 — 2

Рис. 4. Аппроксимация потенциалов взаимодействия в системах И2O•MgO (а) и И2O•ZnO (б) аналитической функцией Букингема 1 - численный эксперимент, 2 - аналитическая функция.

б

Заключение. В настоящей работе вычислены потенциалы взаимодействия молекулы воды с поверхностями кристаллов MgO и ZnO, имеющих гранецентрированную кубическую решетку. Проведено сравнение потенциалов для систем H2O-MgO (рис. 4, а) и H2O-ZnO (рис. 4, б), полученных с помощью квантово-механического метода функционала электронной плотности с гибридным потенциалом B3LYP/6-31G. Определенные потенциалы аппроксимированы функцией Букингема различными подходами: с использованием аналитических вычислений и численных методов в программных пакетах Gnuplot и Origin. Следует отметить, что оба направления имеют свои достоинства и недостатки, но для решения поставленной задачи, а именно нахождения аналитического вида потенциала межмолекулярного взаимодействия для его применения в расчетах методами молекулярной динамики, более всего подходит программный пакет Gnuplot.

В дальнейшем для моделирования взаимодействия большого числа молекул воды над поверхностью кристаллов MgO и ZnO в программе DL POLY [8] планируется взять за основу вид потенциала, полученный в программном пакете Gnuplot.

Литература

1. Егоров Н. В., Денисов В. П., Ермошина М. С. Моделирование взаимодействия заряженной частицы с поверхностью электрода сложной конфигурации // Поверхность. 2005. № 7. C. 60—63.

2. Koch W., Holthausen M. C. A chemist's guide to density functional theory. Ed. 2. Weinheim: Wiley-VCH, 2002. 293 p.

3. Frisch M. J., Trucks G. W., Schlege H. B. e. a. GAUSSIAN 03, Rev. B.05. Pittsburgh, PA: Gaussian, 2003. 596 p.

4. Бедрина М. Е., Егоров Н. В., Клемешев В. А. Моделирование наноструктур на высокопроизводительном вычислительном комплексе // Вестн. С.-Петерб. ун-та. Сер. 4: Физика, химия. 2010. Вып. 4. C. 136-140.

5. Бедрина М. Е., Егоров Н. В., Куранов Д. Ю., Семенов С. Г. Расчет фталоцианинатов цинка на высокопроизводительном вычислительном комплексе // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2011. Вып. 3. С. 13-21.

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

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

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

8. Todorov I. T., Smith V., Trachenko K., Dove M. T. DL POLY 3 new dimensions in molecular dynamics simulations via massive parallelism // J. of materials chemistry. 2006. Vol. 16. P. 1911-1918.

Статья принята к печати 26 апреля 2012 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.