Статья опубликована в №88 (декабрь) 2020
Разделы: Физика, Химия
Размещена 16.12.2020. Последняя правка: 19.12.2020.
Просмотров - 982

ПРИМЕНЕНИЕ МЕТОДОВ МОЛЕКУЛЯРНО-МЕХАНИЧЕСКОГО РАСЧЕТА ДЛЯ ПОСТРОЕНИЯ И ОПРЕДЕЛЕНИЯ СВОЙСТВ КЛАСТЕРОВ МЕТАНА

Голубев Владимир Константинович

Кандидат физико-математических наук, доцент

Нижний Новгород; Университет Людвига-Максимилиана, Мюнхен

Независимый эксперт; приглашенный ученый

УДК 544.277.6

Введение

Несвязанное дисперсионное взаимодействие молекул метана в его димере изучалось довольно основательно с привлечением различных методов молекулярных и квантово-химических расчетов. Необходимым фактором при проведении такого рода квантово-химических расчетов является обязательный учет корреляционных эффектов. Как известно мономер метана обладает нулевым дипольным моментом, что сводит к минимуму диполь-дипольные и другие эффекты второго порядка, тем самым подчеркивая основной вклад дисперсии в энергию взаимодействия двух и более мономеров. Дисперсионная энергия таких систем уже может быть рассчитана с использованием умеренно больших базисных наборов и включением корреляции, оцененной с помощью метода Меллера-Плессета первого порядка (MP).

В работе [1] для изучения дисперсионного взаимодействия в димере метана использовалось несколько расчетных методов. Прежде всего, проводились расчеты с помощью метода Меллера-Плессета второго порядка (MP2) с использованием базисного набора 6-311G(2d,2p). Основными результатами этих расчетов стали энергия взаимодействия составляющих димер мономеров и оптимальное расстояние между ними. Выполненные таким образом расчеты являлись наиболее точными и в последующем они служили эталоном сравнения для всех других используемых методов расчета. К этим другим методам относились полуэмпирические методы MNDO, AM1 и PM3 и молекулярно-механические силовые поля MM2, MM3 и AMBER. В результате проведенных расчетов было получены значения расстояний между атомами углерода r₀, соответствующие наиболее значительной энергии взаимодействия E_i для расположения мономеров голова к голове. Полученные значения r₀ и E_i составили 3.7 Å и 0.42 ккал/моль для уровня теории MP2/6-311G(2d,2p) и 3.8 Å и 0.30 ккал/моль для того же метода с учетом оценки величины ошибки суперпозиции базисного набора (BSSE). Для молекулярно-механических методов MM2, MM3 и AMBER значения r₀ и E_i составили соответственно 3.5 Å и 0.91 ккал/моль, 3.7 Å и 0.56 ккал/моль и 3.7 Å и 0.46 ккал/моль.

В работе [2] потенциалы межмолекулярного взаимодействия димеров метана рассчитывались методами Хартри-Фока, Меллера-Плессета, связанных кластеров и теории функционала плотности. Сравнительные расчеты выполнялись в основном с использованием димера метана с симметрией D_3d и ориентацией мономеров голова к голове. Основные, наиболее показательные результаты выполненных расчетов показаны на рис. 1. Из результатов, приведенных на рисунке, становится понятно, что при использовании высокоточного метода CCSD(T) точность результатов существенно зависит от размеров используемого базисного набора. Такие же методы, как метод Хартри-Фока и метод теории функционала плотности вообще оказываются не в состоянии описать дисперсионное взаимодействие мономеров метана. Наиболее точный результат для энергии взаимодействия E_i = -0.495 ккал/моль получен для расстояния r₀ = 3.5 Å с использованием  уровня теории CCSD(T)/aug(df,pd)-6-311G**.

 

Рис. 1. Расчетные потенциалы межмолекулярного взаимодействия в димере метана D_3d, полученные в работе [2] с использованием нескольких базисных наборов в методе CCSD(T) (слева) и полученные с использованием нескольких различных методов для базисного набора aug(df, pd)-6-311G** (справа).

В работе [3] проводились расчеты для 12 симметричных ориентаций димера метана. Эти ориентации показаны на рис. 2. Расчеты на уровне теории MP3/aug(df,pd)-6-311G** проводились с учетом оценки величины ошибки суперпозиции базисного набора. В результате был получен модельный потенциал exp-6 типа, учитывающий анизотропию H/H взаимодействий.

В работе [4] расчеты для 11 из приведенных на рис. 2 12 ориентаций димера метана проводились на уровне теории MP2/6-311G(2df,2pd). Для проверки возможности увеличения точности без существенного увеличения ресурсоемкости проводились также более высокоуровневые расчеты с использованием MP4/6-311G(2df,2pd) и MP4/aug-cc-pVTZ. На основании результатов выполненных расчетов был получен модельный потенциал в форме Морзе. Дополнительно были проведены расчеты для трех конфигураций тримера метана, двух ориентаций тетрамера и одной ориентации пентамера.

 

Рис. 2. Двенадцать ориентаций димера, используемые в работе [3] при построении модельного потенциала для метана.

В работе [5] выполнялись расчеты димера метана на еще более высоком, практически предельном, уровне теории. Функция энергии взаимодействия дисперсионно-связанного димера метана в конформации D_3d с минимальной энергией рассчитывалась с использованием метода связанного кластера с одиночным, двойным и пертурбативным тройным возбуждением [CCSD(T)] в пределе полного базисного набора (CBS). Потенциальная кривая CCSD(T)/CBS была рассчитана с использованием корреляционно-согласованных базисных наборов поляризованной валентности (cc-pVXZ и aug-cc-pVXZ, X = D, T, Q) и приведена к полной потенциальной кривой с помощью нескольких аналитических экстраполяционных методов. Таким образом, предельно точный результат для дисперсионной энергии взаимодействия мономеров в димере метана E_i = -0.510 ккал/моль был получен для расстояния r₀ = 3.63 Å.

В работе [6] проводилось построение и изучение кластеров метана размерами до сорока молекул (CH₄)_n (n ⩽ 40). Использовалось два потенциала: один - хорошо известный потенциал OPLS, а другой - потенциал Морзе, полученный на основе квантово-химических расчетов димера метана. Поиск глобальных минимумов кластеров метана проводился эвристическим методом в сочетании с геометрическими возмущениями. Последовательности структурного роста кластеров на основе потенциала OPLS и потенциала Морзе были исследованы с константами их вращения и структурными ядрами, образованными внутренними молекулами. Геометрии глобального минимума большинства кластеров на основе потенциала OPLS с n ⩽ 30 аналогичны геометриям соответствующих кластеров, основанных на потенциале Морзе. Для всех кластеров получены их наиболее низкие потенциальные энергии и в результате показано, что взаимодействие молекул в кластерах на основе потенциала OPLS является более интенсивным по сравнению с кластерами на основе потенциале Морзе.

В арсенал популярной, используемой в основном в учебных целях, квантово-химической программы HyperChem [7] входят 4 молекулярно-механических потенциала, или силовых поля, такие как MM+, AMBER, OPLS и BIO+. Силовое поле ММ+ разрабатывалось для широкого круга органических молекул и является базовым силовым полем молекулярной механики в программе HyperChem. Прототипом его является силовое поле ММ2, которое было существенно модифицировано разработчиками программы HyperChem путем уточнения вида потенциальных функций и добавления параметров для новых типов атомов. Можно специально отметить, что энергия ван-дер-ваальсовских взаимодействий вычисляется в этом силовом поле с использованием потенциала "6-ехр".

Силовое поле AMBER разрабатывалось, прежде всего, для расчета белков и нуклеиновых кислот. Потенциальные функции энергии деформации связей и валентных углов ограничиваются в поле AMBER только гармоническим приближением. Для ван-дер-ваальсовских взаимодействий в этом силовом поле используется потенциал "6-12".

Силовое поле OPLS также было предназначено для моделирования белков и нуклеиновых кислот. Однако его особенность состоит в том, что параметры невалентных взаимодействий были получены на основании тщательного моделирования соответствующих низкомолекулярных соединений. Аналитический вид потенциальных функций аналогичен таковому для поля AMBER, за исключением того, что параметры ван-дер-ваальсовскихих взаимодействий вычисляются по несколько иным формулам.

Силовое поле BIO+ является развитием силового поля СНАRМ, созданного для моделирования макромолекул. Аналитический вид потенциальных функций этого силового поля аналогичен таковому для поля AMBER. Некоторые особенности, которые могут быть отнесены к недостаткам, заключаются в том, что в этом поле не производятся расчеты взаимодействия неподеленных электронных пар и водородного связывания.

В данной работе поставлена задача рассмотреть возможность построения и определения свойств не только малых, но и средних по классификации [8] молекулярных кластеров метана с использованием силовых полей молекулярной механики, задействованных в программе HyperChem. При этом на этом начальном этапе исследования несколько большее внимание уделяется более популярному для решения задач расчета структуры и взаимодействия небольших молекул силовому полю ММ+. Также на этом этапе исследования для всех силовых полей  используются настройки, принятые в программе по умолчанию, и никаких экспериментов с изменением настроек программы и изменением параметров используемых силовых полей не проводится.

Результаты и обсуждение

Сначала рассмотрим результаты расчетов различных кластеров метана, полученные с использованием силового поля ММ+. Для димеров в этом случае были получены две стационарные конфигурации, которые соответствуют приведенным на рис 2 ориентациям H и G. Эти димеры показаны на рис. 3, а рассчитанные для них значения энергии E_t и расстояния между атомами углерода d_t составляют соответственно -0.9126 ккал/моль и 3.5019 Å (H) и -0.8790 ккал/моль и 3.5371 Å (G). Отметим также, что для мономера, то есть изолированной молекулы метана, значение E_t было нулевым.

 

Рис. 3. Стационарные конфигурации димеров метана при их расчете с использованием силового поля ММ+: сверху – ориентация H и снизу – ориентация G из рис. 2.

Форма расчетного потенциала межмолекулярного взаимодействия в наиболее стабильном димере показана на рис. 4. В данном случае геометрия молекул, оптимизированная для стационарного состояния димера, оставалась неизменяемой при расчетах с различными значениями расстояния между атомами углерода d.

 

Рис. 4. Расчетный потенциал межмолекулярного взаимодействия в димере метана с ориентацией H при использовании силового поля ММ+.

Последующие кластеры, от тримера до гексамера, имели определенную стабильную геометрию их углеродных остовов и показаны на рис. 5 и 6. Формы тримера и тетрамера были близки к равностороннему треугольнику и правильному тетраэдру (рис. 5). Расстояния между атомами углерода в тримере имели среднее значение 3.6815 Å, незначительно отличаясь от него в четвертом после точки знаке. Длины ребер тетраэдра в тетрамере имели среднее значение 3.8497 Å и незначительные отличия от него уже во втором знаке.

 

Рис. 5 Кластеры, состоящие из трех (слева) и четырех (справа) молекул метана.

Формы пентамера и гексамера напоминали тригональную бипирамиду и октаэдр (рис. 6). Длины ребер тетраэдра в бипирамиде имели среднее значение 3.8840 Å, а в октаэдре – 3.9224 Å.

 

Рис. 6. Кластеры, состоящие из пяти (слева) и шести (справа) молекул метана.

Дальнейшее построение более крупных кластеров осуществлялось путем добавления к уже построенным экземплярам дополнительных молекул метана и оптимизации образующихся структур. Для улучшения этого процесса использовался своего рода механизм утряски, когда на короткое время запускался процесс оптимизации с использованием какого-либо другого силового поля. За окончательный вариант бралась рассчитанная структура с минимальной энергией.

По аналогии с подходом, предложенным в работе [9], результаты выполненных для всех кластеров расчетов представлялись в виде удельной (на одну молекулу) энергии полного разложения кластера из n молекул E_d(n) = E_t(1) - E_t(n)/n, то есть значение E_d для кластера из n молекул определялось как разность между энергией молекулы (иными словами, кластера из одной молекулы) и энергией кластера из n молекул, деленной на n. Результаты выполненных таким образом расчетов представлены на рис. 7, а вид кластера максимального размера показан на рис. 8.

 

Рис. 7. Влияние числа молекул на удельную энергию разложения кластера метана при использовании силового поля ММ+.

 

Рис. 8. Кластер максимального размера из 800 молекул при расчете с использованием силового поля ММ+.

С использованием реализованного в программе метода QSAR (Количественное соотношение структура-свойство) были выполнены расчеты объема и площади внешней поверхности построенных кластеров. Эти результаты для кластеров размером до 400 молекул представлены на рис. 9, 10.

 

Рис. 9. Влияние числа молекул на объем кластера метана.

 

Рис. 10. Влияние числа молекул на площадь внешней поверхности кластера метана.

В предположении, что форма кластеров является сферической, на основании результатов для объема и площади были определены эффективные радиусы кластеров, а зависимости их от числа молекул показаны на рис. 11. С учетом слабого предположения о сферичности кластеров и ограниченной точности метода QSAR полученные результаты по объему и плотности следует признать достаточно приемлемыми оценками этих величин.

 

Рис. 11. Влияние числа молекул на эффективные радиусы кластера метана, определенные на основании его объема (ромбы) и площади внешней поверхности (квадраты).

С использованием полученных результатов по объему и известной массы кластеров были также сделаны оценочные расчеты плотности материала кластера в зависимости от его размера. Эти результаты представлены на рис. 12. Можно отметить, что согласно данным, приведенным в работе [3], плотность жидкого метана составляет 0.424 г/см³. В жидком состоянии метан пребывает в диапазоне температур 91-111 К [10], а расчетные результаты можно отнести к температуре 0 К.

 

Рис. 12. Влияние числа молекул на плотность метана в кластере.

Аналогичные результаты по расчету различных кластеров метана были получены с использованием силового поля OPLS. Для димеров в этом случае были получены четыре стационарные конфигурации, которые в порядке увеличения их энергии соответствуют приведенным на рис 2 ориентациям L, G, H, K. Рассчитанные для этих димеров значения энергии E_t и расстояния между атомами углерода d_t составляют соответственно -0.3763 ккал/моль и 3.8914 Å (L), -0.3756 ккал/моль и 3.7779 Å (G), -0.3752 ккал/моль и 3.7779 Å (H) и -0.3644 ккал/моль и 3.9417 Å (K). Отметим также, что для мономера, то есть изолированной молекулы метана значение E_t составляло 0.00022 кДж/моль. Результаты выполненных для всех кластеров расчетов представлены на рис. 13, а вид кластера максимального размера из 180 молекул показан на рис. 14.

 

Рис. 13. Влияние числа молекул на удельную энергию разложения кластера метана при использовании силового поля OPLS.

 

Рис. 14. Кластер из 180 молекул при расчете с использованием силового поля OPLS.

Силовое поле OPLS использовалось в работе [6] для построения и изучения кластеров, содержащих до 40 молекул метана. Полученные в этой работе данные приведены к зависимости E_d(n) и показаны на рис. 15, 16 вместе с результатами расчетов, выполненных в данной работе. По результатам, полученным для малых кластеров с однозначной структурой, видно, что различие в полученных зависимостях обусловлено не различием в компоновке молекул в кластерах, а в различных настройках программ и в различном выборе параметров используемого силового поля OPLS.

 

Рис. 15. Влияние числа молекул на удельную энергию разложения кластера метана при использовании силового поля OPLS в данной работе (ромбы) и в работе [6] (квадраты).

 

Рис. 16. Отношение значений E_d, полученных при использовании силового поля OPLS в данной работе (1) и в работе [6] (2).

Результаты по расчету различных кластеров метана до максимального размера n = 180, выполненные с использованием силового поля AMBER показаны на рис. 17. Для димеров в этом случае были получены две стационарные конфигурации, которые соответствуют приведенным на рис 2 ориентациям H и G. Рассчитанные для этих димеров значения энергии E_t и расстояния между атомами углерода d_t составляют соответственно -0.4971 ккал/моль и 3.6321 Å (H) и -0.4897 ккал/моль и 3.7779 Å (G). Отметим также, что для мономера, то есть изолированной молекулы метана значение E_t составляло 0.00022 кДж/моль.

 

Рис. 17. Влияние числа молекул на удельную энергию разложения кластера метана при использовании силового поля AMBER.

Результаты по расчету различных кластеров метана до максимального размера n = 180, выполненные с использованием силового поля BIO+ показаны на рис. 18. Для димеров в этом случае были получены четыре стационарные конфигурации, которые в порядке увеличения их энергии соответствуют приведенным на рис 2 ориентациям H, G, L, K. Рассчитанные для этих димеров значения энергии E_t и расстояния между атомами углерода d_t составляют соответственно -0.2280 ккал/моль и 3.7671 Å (H), -0.2279 ккал/моль и 3.7711 Å (G), -0.2235 ккал/моль и 3.8770 Å (L) и -0.2136 ккал/моль и 3.9139 Å (K). Отметим также, что для мономера, то есть изолированной молекулы метана, значение E_t составляло 0.0790 ккал/моль.

 

Рис. 18. Влияние числа молекул на удельную энергию разложения кластера метана при использовании силового поля BIO+.

Из сопоставления полученных расчетных результатов для разных силовых полей с результатами точного расчета дисперсионной энергии взаимодействия мономеров в димере метана E_i = -0.510 ккал/моль и межатомного расстояния r₀ = 3.63 Å [5] можно заключить, что наиболее точное согласие дает результат, полученный с использованием силового поля AMBER. Полученные для этого поля результаты дают значения E_i = -0.4972 ккал/моль и r₀ = 3.6321 Å.

Для энергии отрыва одной молекулы от кластера из n молекул получаем выражение E_s(n) = n·[E_d(n)-E_d(n-1)] + E_d(n-1), из которого следует, что при значительном увеличении размера кластера можно оценить значение эффективной энергии (теплоты) сублимации вещества, составляющего кластер, как E_s = E_d. Попробуем сделать такую оценку на основании полученных расчетных результатов и сопоставить результат с известными экспериментальными данными для энтальпии сублимации твердого метана, приведенными в работе [10].

Максимальное число молекул в кластере при расчете с использованием силового поля MM+ составляло 800 (рис. 7), тогда как при расчете с использованием силового поля AMBER только 180 (рис. 17). Используя пропорцию для результатов рис. 7, получим для значения E_d, соответствующего 800 молекулам на рис. 17, ориентировочную оценку 9.51 кДж/моль. В связи с тем, что размеры кластера являются еще недостаточно большими для осуществления предельного перехода к твердому веществу, да и структура кластера является не кристаллической, а аморфной, данное значение можно только рассматривать как вполне корректную нижнюю оценку для энергии сублимации твердого метана при нулевой температуре.

 

Рис. 19. Результаты измерения теплоемкости метана в зависимости от температуры, полученные в работе [11].

Можно также оценить, насколько близко эта оценка будет соответствовать реальным экспериментальным значениям энтальпии сублимации твердого метана, полученным в температурном диапазоне 50-90 K [10]. Эти значения, приведенные по результатам пяти групп исследователей, лежат в пределах 9.2-10.0 кДж/моль. Используем для оценки приведенную на рис. 19 зависимость теплоемкости твердого метана от температуры из работы [11]. Наблюдаемые на этой зависимости особенности обусловлены фазовым превращением в твердом состоянии при температуре около 20 K и плавлением при температуре около 90 K.

Таким образом, с учетом увеличения энтальпии метана при увеличении температуры полученная расчетная нижняя оценка энтальпии твердого метана смещается от 9.5 кДж/моль при 0 К до 9.2 кДж/моль при 70 К, что следует признать достаточно хорошим оценочным результатом.

Заключение

На основании полученных в работе результатов рассмотрения возможностей четырех силовых полей из программы Hyperchem для построения и расчетного анализа кластеров метана можно сделать вывод, что все рассмотренные силовые поля в той или иной степени пригодны для построения кластеров и расчетов дисперсионного взаимодействия неполярных молекул.

Для всех силовых полей построены зависимости энергии разложения кластеров в зависимости от их размеров. Получено, что для силового поля ММ+ значения этой энергия несколько выше, чем для остальных силовых полей.

В случае более подробного рассмотрения результатов расчетов для димера метана можно заключить, что результаты, полученные с использованием силового поля AMBER, наиболее близко согласуются с результатами высокоточных квантово-химических расчетов.

Выполненные оценки плотности кластеров различного размера дают ее зависимость от числа молекул в кластере и вполне реалистичные предельные значения при увеличении числа молекул.

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

1. Metzger T.G., Ferguson D.M., Glauser W.A. A computational analysis of interaction energies in methane and neopentane dimer systems // J. Comput. Chem. - 1997. - Vol. 18, No. 1. - P 70-79.
2. Tsuzuki S., Uchimaru T., Tanabe K. Intermolecular interaction potentials of methane and ethylene dimers calculated with the Møller-Plesset, coupled cluster and density functional methods // Chem. Phys. Lett. - 1998. - Vol. 287, Iss. 1-2. - P. 202-208.
3. Tsuzuki S., Uchimaru T., Tanabe K. A new ab initio based model potential for methane // Chem. Phys. Lett. - 1998. - Vol. 287, Iss. 3-4. - P. 327-332.
4. Rowley R.L., Pakkanen T. Determination of a methane intermolecular potential model for use in molecular simulations from ab initio calculations // J. Chem. Phys. - 1999. - Vol. 110, No. 7. - P. 3368-3377.
5. Li A. H.-T., Chao S. D. Interaction energies of dispersion-bound methane dimer from coupled cluster method at complete basis set limit // J. Mol. Str. (Theochem). - 2009. - Vol. 897, Iss. 1-3. - P. 90-94.
6. H. Takeuchi. The structural investigation on small methane clusters described by two different potentials // Comput. Theor. Chem. - 2012. - Vol. 986. - P. 48-56.
7. HyperChem Release 8. Tools for Molecular Modeling. - Hypercube, Gainesville, FL, 2007. - 2220 p.
8. Johnston R.L. Atomic and Molecular Clusters. - London: Taylor & Francis, 2002. - 250 p.
9. Golubev V.K. Molecular properties and intermolecular interaction for 1,1-diamino-2,2-dinitroethene and its three cyclic homologs // Procidings of XXI Int. Seminar "New Trends in Research of Energetic Materials". - Pardubice, Czech Republic, 2018. - P. 600-611.
10. NIST Chemistry WebBook: NIST Standard Reference Database Number 69. - Gaithersburg: National Institute of Standards and Technology, 2020. URL: https://webbook.nist.gov/chemistry/
11. Colwell J.H., Gill E.K., Morrison J.A. Thermodynamic properties of CH4 and CD4. Interpretation of the properties of the solids. // J. Chem. Phys. - 1963. - Vol. 39, No. 3. - P. 635-653.

Рецензии:

18.12.2020, 1:19 Мирмович Эдуард Григорьевич
Рецензия: Казалось бы тривиальным - предсказать, что с увеличением числа молекул в кластере интегральные параметры должный выйти на логистическую кривую, однако последовательность этого явления, испытание для этой цели специализированного программного комплекса моделирования, изящная реализация расчётных экспериментов, логическая "смесь" с вплетением обзорной информации, "перекличка" в ней недавней авторской работы [9] с ранними работаvи [8] и [11] обеспечивают безусловную научную основательность как постановки задачи по использованию такого инструментария к моделированию кластерной химии, так и полученные результаты, продемонстрированные на высоко профессиональном уровне. Хотелось бы чуть больше интерпретации к некоторым замечательным иллюстрациям, например, о потенциальной яме на рис. 4. Придирчивый рецензент не нашёл и в оформлении статьи недостатков. Вывод. Данная работа - одно из украшений журнала, и рекомендуется к печати.

18.12.2020, 11:20 Ощепкова Юлия Игоревна
Рецензия: Статья выполнена на высоком уровне, написана хорошим языком. Рекомендую статью к публикации.



Комментарии пользователей:

Оставить комментарий