CC BY
116
УДК 539.194
Т. А. Андреева, М. Е. Бедрина
Вестник СПбГУ. Сер. 10. 2015. Вып. 1
ВЛИЯНИЕ ГИБРИДНЫХ ПОТЕНЦИАЛОВ МЕТОДА DFT НА РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ ЖИДКОКРИСТАЛЛИЧЕСКОЙ ФАЗЫ ВЕЩЕСТВА*)
Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7/9
Квантово-механическим методом функционала электронной плотности (DFT) изучены геометрия, электронная структура димеров и тримеров молекулы 4-га-пентил-4'-цианобифенил (CB5) в программном пакете Gaussian 09. Исследовано влияние гибридных потенциалов различного вида на результаты расчетов. Показано, что процесс сходимости с оптимизацией геометрии надежнее реализуется при использовании потенциала B3LYP и затраты процессорного времени при этом меньше. Энергии связи в ассоциатах, полученные разными методами DFT, примерно одинаковы. Библиогр. 12 назв. Ил. 2. Табл. 3. Ключевые слова: квантовая химия, B3LYP, DFT, M06, жидкие кристаллы.
Т. A. Andreeva, M. E. Bedrina
THE INFLUENCE OF HYBRID POTENTIALS OF THE DFT METHOD ON THE RESULTS OF STUDY OF LIQUID CRYSTAL PHASE OF A SUBSTANCE
St. Petersburg State University, 7/9, Universitetskaya embankment, St. Petersburg, 199034, Russian Federation
Geometry and electronic structure of the dimers and trimers in CB5 molecule were studied by quantum-mechanical DFT method in the software package Gaussian 09. The influence of different types hybrid potentials on the results of the calculations was investigated. It is shown that the convergence process in the optimization of the geometry is more reliably implemented using B3LYP method. Binding energy in the associates obtained by different DFT methods are approximately the same. Bibliogr. 12. Il. 2. Tables 3.
Keywords: quantum chemistry, B3LYP, DFT, M06, liquid crystals.
Введение. Теория межмолекулярных взаимодействий в настоящее время - широко развивающаяся область науки, с ней приходится сталкиваться при рассмотрении не только актов соударения инертных атомов, но и процессов взаимодействия сложных биологически активных молекул и макроскопических тел, при нахождении равновесной геометрии кристаллов и во многих других процессах. Расчеты потенциальных кривых даже в случае простейших систем представляют сложную задачу. Слабые ван-дер-ваальсовы взаимодействия очень чувствительны к выбору метода расчета в целом, в частности к выбору потенциала [1].
Выбранная в качестве объекта исследования молекула цианобифенила и ее ассо-циаты являются одними из типичных представителей жидких кристаллов. Квантово-
Андреева Татьяна Анатольевна — аспирант, e-mail: andreeva_t90@mail.ru
Бедрина Марина Евгеньевна — кандидат химических наук, доцент; e-mail: m.bedrina@mail.ru Andreeva Tatiana Anatolyevna — post-graduent student; e-mail: andreeva_t90@mail.ru Bedrina Marina Evgenievna — candidate of chemical sciences, associated professor; e-mail: m.bedrina@mail.ru
*) Работа выполнена при финансовой поддержке Санкт-Петербургского государственного университета (тема № 13-01-00150).
механическое исследование процесса образования жидкокристаллической (ЖК) фазы на атомно-молекулярном уровне с учетом всех типов межмолекулярных сил до сих пор не проводилось. Решение этой задачи требует применения сложных, современных, расчетных методов и больших затрат вычислительного ресурса. В работе были рассмотрены ассоциаты (димеры, тримеры) молекулы 4-п-пентил-4'-цианобифенила (CB5) и найдены комплексы с наиболее энергетически выигрышной структурой [2, 3]. Компьютерное моделирование проводилось с использованием различных гибридных потенциалов метода функционала электронной плотности [4]. Анализ полученных результатов позволяет выявить оптимальные для учета ван-дер-ваальсовых взаимодействий потенциалы метода функционала электронной плотности (Density Functional Theory, DFT) [5-8].
Математическая модель. Математический аппарат опирается на основные методы квантовой механики и методы квантовой химии. В основе метода Хартри-Фока лежит представление о волновой функции как характеристике состояния квантовой системы. Зная волновую функцию, можно дать точный ответ о строении, энергии и химических свойствах исследуемого соединения. Для этого необходимы полный учет энергии коррелированного движения электронов и представление атомных ор-биталей, не содержащее погрешности базисного набора. В настоящий момент такие расчеты возможны только для самых простых молекул. В зависимости от приближений, используемых в методе Хартри-Фока, встречаются различные модификации расчетных программ. Однако метод не позволяет учитывать эффект корреляции электронов для таких сложных систем как жидкие кристаллы, для которых свойственны слабые эффекты взаимодействия между молекулами при образовании ЖК фазы. Альтернативой таким методам служит подход, основанный на применении теории функционала плотности DFT. Полную энергию квантовой системы можно записать в виде
¿total = T + Ene + J + K + Enn.
С учетом приближения Борна-Оппенгеймера энергия межъядерного отталкивания постоянна, Enn = const. Энергия притяжения электронов к ядрам Ene и кулонов-ская энергия отталкивания электронов J могут быть выражены через электронную плотность основного состояния системы p(r):
N M z N Z
E„e(fi(r)) = 4>i{r)-^4>i{r)dr = / -уp(r)dr,
\ Я \Ф«г) Г рт^ \Ь(г)\2 dr1*2 = \ Я
Начало применению БЕТ методов положило внедрение в расчетную схему орби-талей, предложенное Коном и Шамом. Основная идея теории состоит в разделении функционала кинетической энергии Т на две части, первая вычисляется точно с использованием формально построенных орбиталей, отвечающих системе невзаимодействующих электронов Т$, вторая представляет собой поправочный член - коррекцию Тс:
Т (р)= Тв + Тс.
Коррекция к кинетической энергии Тс в методе Кона-Шама просто включается в обменный член К, и, таким образом, полная электронная энергия в теории БЕТ вычисляется следующим образом:
£врт(р) = Ts (р) + Епе(р) + J (р) + Ехс (р),
где Тз - кинетическая энергия, рассчитываемая с помощью детерминанта Слэтера (Б), составленного из орбиталей Кона-Шама; Ехс - корреляционно-обменный функционал, причем обменный член К дает больший вклад. Различные БЕТ методы отличаются друг от друга выбором формы Ехс функционала, а также наличием наборов подгоночных параметров. Обычным подходом является разделение Ехс на обменный Ех и корреляционный Ес потенциалы:
Ехс (р) = Ех(р) +Ес(р) = ! р(т)ех (р(г))йт + ^ р(т)ес(р(т))йт.
В последнее время стали развиваться гибридные методы, в которых в выражение для Ехс вносится вклад, рассчитываемый так же, как обменная энергия в методе Хартри-Фока, например трехпараметровый функционал Бекке (В3):
ЕВ3 = (1 - а)Е^}иА + аЕнр + ЬАЕВ88 + Е+ сАЕ%СА,
здесь хартри-фоковская обменная энергия Ен с вкладом а входит в выражение для обменного функционала. Последний член суммы есть градиентная коррекция функционала Ес . В зависимости от вида поправки получается та или иная конкретная реализация гибридного БЕТ метода. В наиболее используемом в настоящее время БЕТ методе вместо Е^°А применяется корреляционный функционал Ли, Янга и Парра LYP - B3LYP. Так как функционал LYP содержит и локальный, и нелокальный члены, то корреляционно-обменный функционал метода B3LYP имеет вид
Е^Р = (1 - а)ЕхзОА + аЕнр + ЬАЕВ88 + (1 - с)Е+ сАЕ^Р,
где а, Ь и с - константы, которые были подобраны Бекке с использованием надежных экспериментальных данных для репрезентативного набора сравнительно простых химических соединений.
В 2008 г. в работе Труяра и Зао [9] был предложен гибридный потенциал, специально калиброванный таким образом, чтобы хорошо описывать слабые ван-дер-ваальсовы взаимодействия, чувствительные к выбору метода расчета в целом, в частности к выбору потенциала. Локальные части функционала М06 зависят от трех переменных: спиновой плотности ра, градиента снижения спиновой плотности ха и спиновой плотности кинетической энергии та:
, , оссир
У/У о 1 ,™у ,2
Ха = 4/3 , <Т = а,[3, Та = - ■
ра г
Функционал М06 связывает с рабочей переменной ха две рабочие функции 7 и х:
2та „ „ 3...
■7Г
" -Ср, Ср = —(67Г2)2/3,
р5/з " * 5
ра
Ло , ¿1х2а + ¿2га ¿зхЬ + ¿4x1 га +
НХа, ¿<т) = -Г Н--07-7- +
1(Ха ) 72 (Ха ,га) 73 (Ха ,га)
1 (ха, ¿а) = 1 + а(х1 + га), а обозначает компоненту вдоль произвольного пространства, зафиксированного осью
2
электронного углового момента вращения. Запишем гибридный мета-функционал М06 так:
ЕМ06 = £ /(^(Ра, VЮ)+£х.ОЛЬх(Ха
а
здесь ЕрВЕ(ра, Vра) - плотность обменной энергии модели РВЕ, е!^^ - локальное приближение спиновой плотности
/ 3 \ 1/3 =-3/2 (¿)
f (wa) - фактор повышения спиновой кинетической плотности энергии
m
f ) =
а
т
(Ыа ) = У^ а^юа,
г=0
где переменная юа - функция 1а, а 1а - функция спиновой плотности кинетической энергии та и спиновой плотности ра:
+ _ 1 т ЬЗБЛ
— + — а
■и>а . , 1 ' <т I
tа + 1 Та
ЬЗБА _ ^ ^4/3
Для компьютерного моделирования ЖК фазы были выбраны следующие гибридные потенциалы метода DFT - метод B3LYP и модификации метода M06: M06, М062Х, M06L. Расчеты проводились в программном пакете Gaussian 09 [10] с использованием высокопроизводительного вычислительного комплекса факультета прикладной математики-процессов управления СПбГУ. Для решения поставленной задачи во всех применяемых методах DFT необходимо найти корреляционно-обменную энергию. Проблема состоит в том, что вид корреляционно-обменного функционала неизвестен, а потому необходимо искать его приближенное описание и вычислять численно. Для этого выбирается пространственная сетка интегрирования. Сетка для Ехс играет такую же роль, как базис для остальных вкладов энергии в Etotal.
Результаты. Были построены физические модели структуры ассоциатов молекулы CB5, которые могут образовываться в ЖК фазе. Минимизируя функционал полной энергии, проведена оптимизация геометрии всех построенных модельных структур. Наиболее энергетически выгодные структуры в дальнейшем были рассчитаны с помощью различных потенциалов метода DFT и приведены в табл. 1. Все расчеты проводились с использованием валентно-расщепленного базиса 6-31G. При изучении процесса образования ассоциатов рассматриваются разные варианты задания начального положения молекул. Сложность нахождения равновесной геометрии, которая вычисляется методом наискорейшего спуска, заключается в существовании нескольких локальных минимумов энергии на потенциальной поверхности [11, 12]. Глобальному минимуму отвечает наименьшее значение полной энергии системы. Энергия связи отдельных молекул в комплексе определялась по формуле
ДЕ = E - nEM,
где Е - полная энергия комплекса; Ем - энергия изолированной молекулы; n - число молекул ассоциата.
ю о
Таблица 1. Значения полных энергий Е (в а.е.) и энергий связи АЕ (в ккал/моль) в ассоциатах
Название комплекса ЁВЗЬУР АЕВЗЬУР ЁМ06 АЕмоб -^М062Х ДЁМ062Х ЁМОбЬ АЕмобЬ
-1503.8928 7.83 -1502.7461 8.75 —1503.2433 7.5 -1503.7104 7.07
Б1т2 -1503.8837 7.82 -1502.7462 8.82 —1503.2360 2.87 -1503.7061 4.37
Б1тЗ -1503.8884 5.06 -1502.7489 10.54 -1503.2483 10.6 -1503.7144 9.58
Б1т4 -1503.88372 2.15 -1502.7365 2.74 -1503.2387 4.55 -1503.7061 4.36
В1т5 -1503.8885 5.14 -1502.7459 8.05 -1503.2509 12.25 -1503.7104 7.09
Тпт1 -2255.8396 12.05 -2254.1249 15.79 -2254.8647 11.02 -2255.5842 22.27
Тпт2 -2255.8389 11.59 -2254.1249 15.79 -2254.8722 15.69 -2255.579 19.33
ТптЗ -2255.8368 10.26 -2254.1163 10.38 -2254.8602 8.17 -2255.5711 14.06
Тпт4 -2255.8295 5.7 -2254.1057 3.72 -2254.8604 8.3 -2255.5561 4.58
Тпт5 -2255.8304 6.26 -2254.1171 10.86 -2254.8841 23.18 -2255.5812 20.34
Таблица 2. Минимальные расстояния К (в А) между молекулами в ассоциатах
Название комплекса Метод
ВЗЬУР М06 М062Х М06Ь
В1тег1 2.35 2.28 2.30 2.20
В1тегШ 2.56 2.49 2.55 2.53
БтегГУ 3.21 2.99 2.94 2.94
БтегУ 2.49 2.79 2.79 2.28
Тптег1 2.46 2.24 2.29 2.30
ТптегП 2.33 2.24 2.28 2.37
ТптегШ 2.36 2.26 2.42 2.27
ТптегУ 2.64 2.52 2.38 2.30
Таблица 3. Частоты нормальных колебаний ассоциатов (в см-1)
Название комплекса, частота Метод
ВЗЬУР М06 М062Х М06Ь
Б1тег1, С7-Ш 2295 2269 2376 2269
Б1тег1, С25-К2 2293 2272 2377 2272
БтегШ, С7-Ш 2297 2298 2365 2263
БтегШ, С25-К2 2293 2296 2361 2261
Б1тег1У, С7-Ш 2296 2305 2368 2265
Б1тег1У, С25-К2 2294 2301 2365 2269
БтегУ, С7-Ш 2295 2298 2359 2272
БтегУ, С25-К2 2293 2297 2364 2270
ТптегП, С7-Ш 2295 2296 2360 2266
ТптегП, С25-К2 2294 2305 2370 2265
ТптегП, С43-Ш 2290 2303 2367 2260
ТптегШ, С7-Ш 2294 2299 2361 2266
ТптегШ, С25-К2 2297 2303 2369 2260
ТптегШ, С43-Ш 2296 2305 2366 2263
ТптегУ, С7-Ш 2296 2302 2369 2272
ТптегУ, С25-К2 2296 2301 2365 2261
ТптегУ, С43-Ш 2295 2295 2361 2262
В целом применение функционалов B3LYP и М06 дает одинаковую картину образования димеров и тримеров, т. е. наиболее выгодными энергетически были похожие конфигурации, которые незначительно различались относительным расположением молекул. Но структура ван-дер-ваальсовых комплексов все же лучше передается методами М06, а именно методами M06L, M062X. Метод B3LYP сразу дает более низкие величины энергии (уже на 20-м шаге), но при достижении сходимости (50-й шаг и далее), значения энергии сближаются. Характерный вид кривой сходимости димеров на примере Dimerl представлены на рис. 1, а сама структура комплекса - на рис. 2.
Энергия,
ккал/моль —я— М06
18 : 16 J 1412-io-8 -6: 420:
—2 -I—i-1-1-1-1-1-1->—|—■—i—■—|—■—|-1
0 5 10 15 20 25 30 35
Шаг
Рис. 1. Кривая сходимости для структуры DimerI
С точки зрения сходимости процесса оптимизации геометрии, метод B3LYP оказывается более надежным.
Некоторые различия в геометрии ассоциатов характеризуются различными минимальными расстояниями между молекулами (табл. 2).
Частоты характеристического валентного колебания CN связи в отдельной молекуле CB5 и ее ассоциатов, рассчитанные разными методами, отличаются по величине (табл. 3). Наиболее близки к экспериментальным значениям частоты, вычисленные методом M06L.
Заключение. Квантово-механическим методом DFT изучены геометрия, электронная структура димеров и тримеров молекулы CB5 в программном пакете Gaussian 09. Исследовано влияние гибридных потенциалов различного вида на результаты расчетов. Показано, что процесс сходимости при оптимизации геометрии надежнее реализуется при использовании B3LYP и затраты процессорного времени при этом меньше. Энергии связи в ассоциатах примерно одинаковы, за исключением тримера 5 (Trim5), который за счет лучшего учета ван-дер-ваальсовых взаимодействий иначе ориентируется в пространстве.
B3LYP
N2
Н1
Н20
N1
*
Рис. 2. Структура комплекса Бш1ег1
Литература
1. Каплан И. Г. Введение в теорию межмолекулярных взаимодействий. М.: Наука, Гл. ред. физ.-мат. лит., 1982. 312 с.
2. Андреева Т. А., Бедрина М. Е. Моделирование механизма взаимодействия молекул в жидкокристаллической фазе // Процессы управления и устойчивость: Труды 43-й междунар. науч. конференции аспирантов и студентов / под ред. А. С. Ерёмина, Н. В. Смирнова. СПб.: Издат. Дом С.-Петерб. гос. ун-та, 2012. С. 483-487.
3. Андреева Т. А. Анализ гибридных потенциалов метода функционала электронной плотности // Процессы управления и устойчивость: Труды 44-й междунар. науч. конференции аспирантов и студентов / под ред. Н. В. Смирнова, Т. Е. Смирновой. СПб.: Издат. Дом С.-Петерб. гос. ун-та, 2013. С. 135-141.
4. Кларк Т. Компьютерная химия: практ. руководство по расчетам структуры и энергии молекулы / пер. с англ. А. А. Коркина; под ред. В. С. Майстрюкова, Ю. Н. Панченко М.: Мир, 1990. 383 с. (Clark T. A Handbook of computational chemistry.)
5. Wolfram Koch, Max C. Holthausen. A chemist's guide to density functional theory. Second Edition. Berlin, Germany: Wiley-VSN, 2001. 306 p.
6. Yan Zhao, Donald G. Truhlar. Density functionals with broad applicability in chemistry // Accounts of chemical research. 2008. Vol. 41, N 2. P. 157-167.
7. Райк А. В., Бедрина М. Е. Моделирование процесса адсорбции воды на поверхности кристаллов // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2010. Вып. 3. C. 13-21.
8. Бедрина М. Е., Егоров Н. В., Куранов Д. Ю., Семенов С. Г. Расчет фталоцианинатов цинка на высокопроизводительном вычислительном комплексе // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2011. Вып. 3. C. 13-21.
9. Yan Zhao, Donald G. Truhlar. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals // Theor. Chem. Account. USA. 2008. P. 215-241.
10. Frisch M. J., Trucks G. W., Schlegel H. B. e. a. GAUSSIAN-09, Rev. C.01. Wallingford CT: Gaussian, Inc., 2010.
11. Бедрина М. Е., Егоров Н. В., Клемешев В. А. Моделирование наноструктур на высокопроизводительном вычислительном комплексе // Вестн. С.-Петерб. ун-та. Сер. 4: Физика, химия. 2010. Вып. 4. C. 136-140.
12. Райк А. В., Егоров Н. В., Бедрина М. Е. Моделирование потенциалов межмолекулярного взаимодействия // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2012. Вып. 3. C. 79-87.
References
1. Kaplan I. G. Vvedenie v teoriju mezhmolekuljarnyh vzaimodejstvij (Introduction to the Theory of Intermolecular Interactions). Moscow: Nauka, 1982, 312 p.
2. Andreeva T. A., Bedrina М. Е. Modelirovanie mehanizma vzaimodejstvija molekul v zhidkokristallicheskoj faze (Modeling the mechanism of interaction molecules in the liquid crystal phase). The XLIII annual international conference on Control Processes and Stability (CPS'13). Ed. by A. S. Eremin, N. V. Smirnov. St. Petersburg: Izdat. Dom St. Petersburg University, 2012, pp. 483-487.
3. Andreeva T. A. Analiz gibridnyh potencialov metoda funkcionala jelektronnoj plotnosti (Analysis of the hybrid potential of the density functional theory). The XLIV annual international conference on Control Processes and ¡Stability (CPS'14). Ed. by N. V. Smirnov, T. E. Smirnova. St. Petersburg: Izdat. Dom St. Petersburg University, 2013, pp. 135-141.
4. Clark T. A. Komp'juternaja himija: prakt. rukovodstvo po raschetam struktury i jenergii molekuly (Handbook of Computational Chemistry: A Practical Guide to Chemical Structure and Energy Calculations Hardcover). Moscow: Mir, 1990, 383 p.
5. Wolfram Koch, Max C. Holthausen. A chemist's guide to density functional theory. Second Edition. Berlin, Germany: Wiley-VSN, 2001, 306 p.
6. Yan Zhao, Donald G. Truhlar. Density functionals with broad applicability in chemistry. Accounts of chemical research, 2008, vol. 41, no. 2, pp. 157-167.
7. Raik A. V., Bedrina М. Е. Modelirovanie processa adsorbcii vody na poverhnosti kristallov (Modeling of water adsorption process on crystal surfaces). Vestn. of St. Petersburg University. Serie 10: Applied mathematics, computer science, control processes, 2011, issue 2, pp. 67-75.
8. Bedrina M. E., Egorov N. V., Kuranov D. Y., Semenov S. G. Raschet ftalocianinatov cinka na vysokoproizvoditel'nom vychislitel'nom komplekse (Calculation phthalocyaninatozinc on the high-efficiency computer complex). Vestn. of St. Petersburg University. Serie 10: Applied mathematics, computer science, control processes, 2011, issue 3, pp. 13—21.
9. Yan Zhao, Donald G. Truhlar. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Account. USA, 2008, pp. 215-241.
10. Frisch M. J., Trucks G. W., Schlegel H. B. e. a. GAUSSIAN-09, Rev. C.01. Wallingford CT: Gaussian, Inc., 2010.
11. Bedrina M. E., Egorov N. V., Klemeshev V. A. Modelirovanie nanostruktur na vysokoproizvoditel'nom vychislitel'nom komplekse (Modeling nanostructures on a high performance computing system). Vestn. of St. Petersburg University. Serie 4: Physics, chemistry, 2010, issue 4, pp. 136-140.
12. Raik A. V., Egorov N. V., Bedrina M. E. Modelirovanie potencialov mezhmolekuljarnogo vzaimodejstvija (Modelling intermolecular interaction potentials). Vestn. of St. Petersburg University. Serie 10: Applied mathematics, computer science, control processes, 2012, issue 3, pp. 79-87.
Статья рекомендована к печати проф. Н. В. Егоровым. Статья поступила в редакцию 13 ноября 2014 г.
