CC BY
62
УДК 540.1:532.7
СТРУКТУРА СВЕРХКРИТИЧЕСКОГО МЕТАНОЛА ПО ДАННЫМ КЛАССИЧЕСКОЙ МОЛЕКУЛЯРНОЙ ДИНАМИКИ И МЕТОДА КАРА-ПАРИНЕЛЛО
Н.А. Абакумова1, Е.Г. Одинцова2, В.Е. Петренко2
Кафедра «Химия», ФГБОУВПО «ТГТУ» (1);
ФГБУН «Институт химии растворов им. Г.А. Крестова РАН», г. Иваново (2);
vep@isc-ras.ru
Ключевые слова и фразы: водородная связь; сверхкритическое состояние; методы молекулярной динамики, Монте-Карло, Кара-Паринелло; функция радиального распределения.
Аннотация: Проведен расчет функций радиального распределения и среднего числа водородных связей на молекулу метанола при нормальных, суб- и сверхкритических условиях в рамках классической молекулярной динамики и молекулярной динамики Кара-Паринелло. Показано, что независимые методы моделирования дают близкие результаты. При этом рассчитанные функции радиального распределения хорошо согласуются с экспериментом только при T = 298 К, P = 0,1 МПа, а при высоких температурах и давлениях наблюдаются значительные расхождения с экспериментальными функциями, известными из литературы. С другой стороны, оба метода моделирования очень близко к эксперименту воспроизводят степень водородного связывания в метаноле и ее изменение в зависимости от параметров состояния.
Введение
Интерес к простым спиртам в сверхкритическом (СК) состоянии, резко возросший в последние два десятилетия, во многом обусловлен их огромной практической значимостью. Наиболее популярным является метанол: в сверхкритиче-ском состоянии он широко используется в роли растворителей, реагентов и катализаторов. Подобный интерес связан, прежде всего, с перспективами его практического применения в качестве альтернативного СК-воде растворителя, который проявляет высокую растворяющую способность по отношению к широкому классу соединений и при этом менее агрессивен к аппаратуре и обладает более низкими критическими параметрами.
Структура метанола в СК-состоянии неоднократно была объектом исследований различных экспериментальных методов: ядерного магнитного резонанса, инфракрасной (ИК) и Раман-спектроскопии. Математическая обработка спектральных данных позволила с высокой степенью достоверности определить ряд структурных характеристик СК-метанола: степень водородного связывания, доли связанных молекул и т. д. Однако единственным прямым источником детальной
информации о молекулярной структуре и микродинамике флюида служат методы численного моделирования. С их помощью можно получить пространственновременные параметры, недоступные из физического эксперимента.
Функции радиального распределения (ФРР) - характеристика, в полной мере отражающая зависимость обобщенной структуры флюида от межмолекулярно-го расстояния. В усредненном виде в них заложена вся информация о траекториях движения каждой отдельной частицы, а также обо всех возможных молекулярных конфигурациях. Поэтому корректное воспроизведение ФРР очень важно при оценке предсказательной способности любого метода моделирования.
Экспериментально найденные ФРР СК-метанола впервые были получены в 2000 г. посредством обработки методом Монте-Карло (МК) парциальных структурных факторов, полученных методом дифракции нейтронов [1]. Эксперимент проводился в трех точках фазовой диаграммы: при T = 475 К, P = 73,7 МПа, р = 0,700 г/см3 (субкритическое состояние), при T = 526 К, P = 117,7 МПа, р = 0,700 г/см3 (высокоплотный сверхкритический флюид) и при T = 526 К, P = 14,3 МПа, р = 0,453 г/см3 (низкоплотный сверхкритический флюид).
Как выяснилось, экспериментальные атом-атомные ФРР [1] существенно отличаются от рассчитанных ранее методом молекулярной динамики (МД) [2, 3]. Позднее исследователи неоднократно обращались к моделированию структуры СК-метанола, но ФРР были рассчитаны только в работе [4], и их отличие от экспериментальных оказалось столь же существенно.
Таким образом, несмотря на значительное количество публикаций, в которых рассматриваются структурные характеристики СК-метанола, и на разносторонность методов исследования, информации о его ФРР очень немного. При этом между экспериментальными и полученными в численном моделировании ФРР существуют расхождения. Данные по ФРР метанола при стандартных условиях (с. у.) T = 298 К, P = 0,1 МПа, напротив, обширны, в них постоянно вносятся уточнения [2, 3, 5], и расчетные ФРР очень близко воспроизводят экспериментальные [6] (рис. 1).
Чтобы проверить, в какой мере численное моделирование способно воспроизвести экспериментальные ФРР метанола при нормальных [6] и сверхкритиче-ских [1] условиях, мы провели собственный независимый расчет с использованием двух методов: классической молекулярной динамики и неэмпирической молекулярной динамики Кара-Паринелло (КПМД) для тех же параметров состояния.
Детали моделирования
Моделирование классической МД выполнено в NVT-ансамбле с использованием программного пакета GROMACS 4.5.4 в полноатомном силовом поле OPLS-AA, в кубической ячейке с периодическими граничными условиями, содержащей 216 молекул метанола. После выхода системы на равновесие продолжительность моделирования составляла 200 пс с шагом интегрирования 0,5 фс.
Метод неэмпирической КПМД был реализован в рамках пакета программ CPMD [11]. Расчет проводился в течение 15 пс в каноническом №УТ-ансамбле. В кубическую ячейку с периодическими граничными условиями были заключены 32 молекулы метанола. Орбитали Кона-Шема раскладывались в базисе плоских волн с энергией обрыва 25 Ry. Для описания электрон-ионных взаимодействий использовался ультрамягкий псевдопотенциал Вандербильта с обменно-корреляционным функционалом BLYP. Фиктивная масса электронов составляла 400 а.е. Шаг интегрирования был принят равным 7 а.е., что соответствует ~ 0,69 фс.
ISSN 0136-5835. Вестник ТГТУ. 2013. Том 19. № 1. Transactions TSTU 123
Рис. 1. Атом-атомные функции радиального распределения метанола #00(/'), ённ(г)> #сс(гХ 8со(!"% Яси(г) (Н - гидроксильный водород) при с. у.:
1 - эксперимент + метод МК (обработка данных дифракции нейтронов методом МК) [6]; 2 - эксперимент + метод МК (обработка данных дифракции рентгеновских лучей методом МК) [7]; 3 - эксперимент + метод МД (обработка данных дифракции нейтронов и дифракции рентгеновских лучей методом МД) [7]; 4 - метод МД (потенциал Н1); 5 - метод МД [2]; 6 - метод МД [3]; 7 - метод КПМД [8]; 8 - метод МД [5]; 9 - метод КПМД [5]; 10 - методы квантовой химии + МД [5]; 11 - метод МД, данные авторов; 12 - метод КПМД, данные авторов
Обсуждение результатов
На рисунках 1-4 сопоставлены ФРР метанола, рассчитанные для н. у. и тех же точек суб- и сверхкритической области фазовой диаграммы, в которых проводился эксперимент [1]. Рисунок 1 наглядно иллюстрирует высокую степень точности, с которой методы численного моделирования (кривые 4-12) воспроизводят экспериментальные (кривые 1-3) ФРР метанола при с. у. Положения основных пиков совпадают, как и их взаимные соотношения по высоте. В абсолютных значениях высот пиков и глубины минимумов наблюдаются расхождения, но они незначительны. Расчетные ФРР отличаются от экспериментальных не более чем сами экспериментальные ФРР (полученные посредством обработки экспериментальных данных разными методами численного моделирования, см. подпись к рис. 1) отличаются друг от друга.
Экспериментальные ФРР [1] (см. рис. 2-4) при переходе в суб- и сверхкрити-ческое состояние изменяются, по сравнению с н. у., причем наиболее выраженные изменения наблюдаются при низкой плотности (см. рис. 3). На всех атом-атомных функциях, за исключением янн(г), намного уменьшаются высота первого максимума и глубина первого минимума, а на gcH(r) первый пик вообще редуцируется в плечо. На gHH(r) первый пик, напротив, увеличивается, а при низкой плотности (см. рис. 3) - более чем в 2 раза, по сравнению с н. у. Также при р = 0,453 г/см3 второй пик на gOH(r) увеличивается настолько, что становится выше первого. Отличие между ФРР при высокой плотности и при н. у. не так значительно, и переход из суб- в сверхкритическое состояние при постоянной р = 0,700 г/см3 очень мало влияет на ФРР (см. рис. 2, 4).
Численное моделирование показывает иную картину. Расчетные ФРР (см. подписи к рис. 1 -4) с ростом температуры и давления меняются иначе, чем экспериментальные ФРР Сопера-Ямагучи (см. рис. 1-4, кривые 1). Если на gOO(r) [1] первый пик не только уменьшается по высоте, по сравнению с н. у., но и расширяется в сторону больших расстояний (второй пик, соответственно, сдвигается вправо), то на gOO(r), полученных в численном моделировании, смещения нет, - меняется только высота пиков. На расчетных ФРР ОН в суб- и сверхкритических состояниях первый пик становится ниже, но второй при этом остается таким же, как и при н. у., а на ФРР ОН [1] одновременно с понижением первого пика происходит увеличение второго. Аналогично ведут себя gCO(r), с той разницей, что второй пик на экспериментальных функциях не только увеличивается, но и сливается при этом с первым, редуцируясь в плечо, - на расчетных этого не наблюдается. Характер изменения расчетных ФРР НН противоположен экспериментальным: первый пик не увеличивается, а, напротив, уменьшается. Общие тенденции в поведении gcc(r) и gCH(r) схожи, различия состоят лишь в том, что пики на расчетных функциях ниже, чем на экспериментальных, а положения их незначительно отличаются.
Итак, между ФРР суб- и сверхкритического метанола, рассчитанными в компьютерном эксперименте и полученными из экспериментальных данных, существует ряд принципиальных различий. Ни метод МД, ни метод КПМД не воспроизводят результатов экспериментов (см. рис. 2-4). С другой стороны, между ФРР, рассчитанными с использованием разных методов моделирования, наблюдаются лишь незначительные количественные расхождения. На рисунке 3, помимо данных собственных расчетов, приведены также ФРР, полученные другими авторами при параметрах состояния, близких к рассматриваемой точке (Т = 526 К; P = 14,3 МПа; р = 0,453 г/см3). Как можно видеть, все ФРР, полученные в численном моделировании, однотипны.
ISSN 0136-5835. Вестник ТГТУ. 2013. Том 19. № 1. Transactions TSTU 125
Рис. 2. Атом-атомиые функции радиального распределения метанола при / = 475 К,Р= 73,7 МПа, р = 0,700 г/см3: 1 - эксперимент + метод МК (обработка данных дифракции нейтронов методом МК) [1]; 2 - метод МД, данные авторов; 3 - метод КПМД, данные авторов
126 ISSN 0136-5835. Вестник ТГТУ. 2013. Том 19. № 1. Transactions TSTU
Рис. 3. Атом-атомиые функции радиального распределения метанола при / = 526 К, Р= 14,3 МПа, р = 0,453 г/см3: 1 - эксперимент + метод МК (обработка данных дифракции нейтронов методом МК) [1]; 2 - метод МД, данные авторов; 3 - метод КПМД, данные авторов; 4 - метод МД [4], Т= 523 К, р = 0,478 г/см3; 5 - метод МД [3], Т= 523 К, р = 0,348 г/см3, Р = 10 МПа; б - метод МД [3], Т= 523 К, р = 0,516 г/см3, Р = 20 МПа
ISSN 0136-5835. Вестник ТГТУ. 2013. Том 19. № 1. Transactions TSTU 127
Рис. 4. Атом-атомиые функции радиального распределения метанола при / = 526 К, Р= 117,7 МПа, р = 0,700 г/см3: 1 - эксперимент + метод МК (обработка данных дифракции нейтронов методом МК) [1]; 2 - метод МД, данные авторов; 3 - метод КПМД, данные авторов
Итак, мы получили подтверждение того, что численное моделирование неспособно воспроизвести ФРР сверхкритического метанола [1] - единственные известные экспериментальные ФРР - даже на качественном уровне. В таблице представлены среднее число инв водородных связей (ВС) на молекулу и степень водородного связывания
пнВ(Т, Р) =-----”нв(Г, Р)------,
«НВ(298 К; 0,1 МПа)
рассчитанные с помощью интегрирования я0Н(г) МД и КПМД на интервале [0; 2,4] А (согласно критерию водородной связи г00 < 3,4 А, г0Н < 2,4 А). Они сопоставлены с инв, полученными Сопером-Ямагучи при интегрировании их ФРР с тем же критерием ВС [1]. Как и следовало ожидать, расхождения оказались значительными.
Но существуют экспериментальные данные о среднем числе ВС на молекулу и степени водородного связывания, полученные не только из ФРР, но и непосредственно из спектральных данных [9]. Поэтому мы дополнительно рассчитали ФРР метанола в трех состояниях (см. таблицу), для которых известна экспериментальная степень водородного связывания [9], вычислили из них инв и пнв и сравнили полученные значения. Как можно видеть, рассчитанные нами инв и Пш очень хорошо согласуются с данными Хоффмана-Конради [9]. Это свидетельствует о корректности получаемых в численном моделировании усредненных характеристик водородного связывания и функций радиального распределения. Значения пнв (метод МД) несколько больше, чем пнв (метод КПМД) в силу того, что первый пик на я0Н(г), получаемых методом МД, систематически выше пика на я0Н(г), получаемых методом КПМД (см. рис. 1 -4). Подобное же явление мы наблюдали при моделировании ФРР воды в широком интервале температур и давлений [10]. Оно обусловлено множеством факторов, и одним из определяющих является размер ансамбля (32 молекулы при использовании метода КПМД и 216 - при МД), с которым непосредственно связана степень учета влияния дальнодействующих корреляций. Величины инв и пнв были рассчитаны методом МД не только путем интегрирования я0Н(г), но также и напрямую - по соответствию критерию ВС, в рамках стандартного программного пакета GROMACS. Полученные значения близки между собой (см. таблицу).
Среднее число ВС на молекулу метанола и степень водородного связывания Пт при суб- и сверхкритических параметрах состояния
T, K P, МПа 3 CL О г/ nHB (КПМД, инт.) Лнв (КПМД, инт.) nHB (МД, инт.) й^ТТ (и nHB (МД, критер.) Лнв (МД, критер.) nHB [1] Пнв [1] 1 1 СО
298 0,1 0,787 1,86 1,00 1,85 1,00 1,85 1,00 1,77; 1,12 [2] 1,00 -
475 73,7 0,700 1,18 0,63 1,22 0,66 1,19 0,64 1,60 0,90 -
526 14,3 0,453 0,79 0,42 0,82 0,44 0,80 0,43 1,00 0,56 0,42*
526 117,7 0,700 1,12 0,60 1,15 0,62 1,11 0,60 1,60 0,90 -
573 15,0 0,210 0,33 0,18 0,35 0,19 0,33 0,18 - - 0,21
573 20,0 0,346 0,48 0,26 0,52 0,28 0,50 0,27 - - 0,29
573 35,0 0,488 0,64 0,34 0,70 0,39 0,67 0,36 - - 0,37
* T = 523 К, P = 15 МПа, р = 0,477 г/см3.
Расхождения между расчетными и экспериментальными ФРР авторы [1] прокомментировали следующим образом: с одной стороны, их причиной является несовершенство методов моделирования, в частности модельного потенциала, с другой, - то, что при обработке экспериментальных парциальных структурных факторов неизбежны систематические ошибки, дающие значительный вклад в ФРР. Возможности методов численного моделирования за 12 лет, прошедших с момента публикации [1], существенно расширились, однако, как мы показали, никакого принципиального усовершенствования в плане моделирования ФРР СК метанола достигнуто не было. Авторы [1] имели возможность сравнить свои ФРР с данными двух работ: [2] (рассмотрена область более низких плотностей) и [3] (при сопоставимых с [3] параметрах состояния). Функции радиального распределения [3], а также более поздние данные [4], приведены на рис. 3 в дополнение к полученным нами функциям. Как можно видеть, все расчетные ФРР однотипны. Относительно же систематических ошибок в процедуре обработки, обуславливающих особое поведение ФРР метанола, нужно отметить, что их вклад становится столь весомым только в суб- и сверхкритических условиях. При нормальных условиях, как показано на рис. 1, различий в поведении экспериментальных и расчетных ФРР не наблюдается.
Данные независимых методов численного моделирования очень близки не только на качественном уровне, но и в количественных оценках структурных характеристик сверхкритического метанола. При этом они хорошо согласуются с экспериментальными данными об изменении степени водородного связывания в СК-метаноле с температурой и давлением [10].
Список литературы
1. Yamaguchi, T. The Structure of Subcritical and Supercritical Methanol by Neutron Diffraction, Empirical Potential Structure Refinement, and Spherical Harmonic Analysis / T. Yamaguchi, C. J. Benmore, A.K. Soper // J. Chem. Phys. - 2000. -Vol. 112, Part 20. - P. 8976-8987.
2. Asahi, N. Nuclear Magnetic Resonance and Molecular Dynamics Study of Methanol Up to the Supercritical Region / N. Asahi, Y. Nakamura // J. Chem. Phys. -1998. - Vol. 109, Part 22. - P. 9879-9887.
3. Chalaris, M. Hydrogen Bonding in Supercritical Methanol. A Molecular Dynamics Investigation / M. Chalaris, J. Samios // J. Phys. Chem. B. - 1999. -Vol. 103, Part 8. - P. 1161-1166.
4. Andanson, J.M. Relation between Hydrogen Bonding and Intramolecular Motions in Liquid and Supercritical Methano / J.M. Andanson, P.A. Bopp, J.C. Soetens // J. Mol. Liq. - 2006. - Vol. 129. - P. 101-107.
5. McGrath, M.J. Liquid Structures of Water, Methanol, and Hydrogen Fluoride at Ambient Conditions from First Principles Molecular Dynamics Simulations with a Dispersion Corrected Density Functional / M.J. McGrath, I.-F.W. Kuo, J.I. Siepmann // Phys. Chem. Chem. Phys. - 2011. - Vol. 13, Part 13. - P. 19943-19950.
6. Yamaguchi, T. The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at - 80 °C and +25 °C / T. Yamaguchi, K. Hidaka, A.K. Soper // Mol. Phys. - 1999. - Vol. 96, Part 8. - P. 1159-1168.
7. Hydrogen Bonding and Molecular Aggregates in Liquid Methanol, Ethanol, and 1-Propanol / A. Vrhovsek [et al.] // J. Phys. Chem. B. - 2011. - Vol. 115, Part 46. -P. 13473-13488.
8. Hydrogen Bond Dynamics in Liquid Methanol / M. Pagliai [et al.] // J. Chem. Phys. - 2003. - Vol. 119, Part 13. - P. 6655-6662.
9. Hoffman, M.M. Are There Hydrogen Bonds шп Supercritical Water / M.M. Hoffman, M.S. Conradi // J. Am. Chem. Soc. - 1997. - Vol. 119. - P. 3811-3817.
10. Гурина, Д.Л. Функции радиального распределения суб- и сверхкритиче-ской воды по результатам неэмпирической молекулярной динамики / Д.Л. Гурина, M^. Антипова, В.Е. Петренко // Журн. физ. химии. - 2011. - Т. 85, вып. 5. -С. 885-901.
11. Петренко, В.Е. Распределения времен жизни водородных связей жидкой воды. Результаты классической молекулярной динамики и динамики ^ра-Парринелло / В.Е. Петренко, Н.А. Абакумова, Д.Л. Гурина // Вестн. Тамб. гос. техн. ун-та. - 2012. - Т. 18, № 2. - С. 411-419.
The Structure of Supercritical Methanol According to Classical Molecular Dynamics and Method of Car-Parrinello
N.A. Abakumova1, E.G. Odintsova2, V.E. Petrenko2
Department “Chemistry ”, TSTU (1);
Institute of Solution Chemistry named after G.A. Krestov Russian Academy of Sciences, Ivanovo (2); vep@isc-ras.ru
Key words and phrases: hydrogen bond; Car-Parrinello; methods of molecular dynamics; Monte Carlo; radial distribution function; supercritical state.
Abstract: The calculation of radial distribution functions and the average number of hydrogen bonds per molecule of methanol under normal, sub-and supercritical conditions in terms of classical molecular dynamics and molecular dynamics of the Car-Parrinello have been made. It is shown that independent modeling techniques give similar results. In this case, the calculated function of the radial distribution is in good agreement with experiment only at T = 298 K, P = 0,1 MPa; while at high temperatures and pressures there observed significant discrepancies with the experimental features known from the literature. On the other hand, both modeling methods reproduce the degree of hydrogen bonding in methanol and the change of the parameters depending on the state very close to the experiment.
Struktur des uberkinetiswchen Methanols nach der klassischen Molekulardynamik und der Methode von Car-Parrinello
Zusammenfassung: Es ist die Berechnung der Funktionen der radialen Verteilung und der mittleren Zahl auf das Molekul des Methanols bei den normalen, sub- und uberkritischen Bedingungen in den Rahmen der klassischen molekularen Dynamik und der molekularen Dynamik von Car-Parrinello durchgefuhrt. Es ist gezeigt, dass die unabhangigen Methoden der Modellierung die naheren Resultate geben. Dabei stehen die berechneten Funktionen der radialen Verteilung mit dem Experiment nur bei T = 298 K, P = 0,1 MPa in Einklang, und bei den hohen T emperaturen und den Drucken warden die wesentlichen Unterschiede mit den aus der
Literatur bekannten experimentellen Funktionen betrachtet. Andererseits reproduzieren die beiden Methoden der Modellierung sehr nah zu dem Experiment den Grad der Wasserstoffbindung im Methanol und seine Veranderung je nach den Parametern des Zustandes.
Structure du methanol ultra-critique d’apres les donnees de la dynamique moleculaire classique et la methode de Car-Parrinello
Resume: Est effectue le calcul des fonctions de la repartition radiale et du nombre moyen des liens de l’hydrogene sur la molecule du methanol lors des conditions normales sous- et ultra-critiques dans le cadre de la dynamique moleculaire classique et celle de Car-Parrinello. Est montre que les methodes independantes du modelage donnent de proches resultats. Avec cela les fonctions calculees de la repartition radiale s’accordent bien avec un experiment seulement lors de T = 298 K, P = 0,1 MPa et lors de hautes temperatures et pressions sont observees les divergences considerables avec les fonctions experimentales connues a partir de la litteraure. De l’autre cote, les deux methodes du modelage sont tres proches a l’experiment et reproduisent le degre du liage hydrogenique dans le methanol et son changement en fonction des parametres de l’etat.
Авторы: Абакумова Нина Алексеевна - кандидат химических наук, доцент кафедры «Химия», ФГБОУ ВПО «ТГТУ»; Одинцова Екатерина Геннадьевна -
аспирант лаборатории «ЯМР-спектрография и численные методы исследования жидких систем»; Петренко Валентина Евгеньевна - кандидат физико-математических наук, старший научный сотрудник лаборатории «ЯМР-спектрография и численные методы исследования жидких систем», ФГБУН «Институт химии растворов им. Г.А. Крестова РАН», г. Иваново.
Рецензент: Килимник Александр Борисович - доктор химических наук, профессор, заведующий кафедрой «Химия», ФГБОУ ВПО «ТГТУ».
