Сер. 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
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 г.