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

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

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

ИНФОРМАТИКА

УДК 539.193

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

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

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

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

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

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

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

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

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

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

© А. В. Райк, Н.В. Егоров, М.Е. Бедрина, 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

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.

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 г.