High-pressure elastic properties of polycyclic aromatic hydrocarbons obtained by first-principles calculations

Full Text

ВВЕДЕНИЕ

Изучение углеводородных соединений при высоком давлении является важной задачей геохимии и химии твердого тела. Знание термодинамики полициклических ароматических углеводородов (ПАУ) необходимо для моделирования твердых и жидких природных углеводородных систем, а также систем с участием С-О-Н и более сложных флюидов (Зубков, 2001; Литасов, 2011). ПАУ, возможно, стабильны при высоких давлениях (Чекалюк, 1967; Зубков, 2009; Чанышев и др., 2014b), так как были найдены в составе включений в алмазах и его минералах-спутниках (Кулакова и др., 1982; Каминский и др., 1985; Гаранин и др., 2011), метеоритном веществе (Oro et al., 1971; Becker et al., 1997), а также являются неотъемлемой частью продуктов ударного взаимодействия в космическом пространстве (Allamandola et al., 1989; Tielens, 2008). По данным изучения инфракрасных спектров межзвездного вещества, более 20% углерода во Вселенной может быть представлено ПАУ (Ehrenfreund, Charnley, 2000; Tielens, 2008).

В рамках данной работы рассмотрены структуры бензола (стабильные при низких температурах и высоком давлении), а также несколько кристаллических структур ПАУ от нафталина до коронена, стабильных при стандартных условиях (0.1 МПа и 298.15 K). При этом учитывались различные конфигурации веществ в зависимости от соединения бензольных колец. Предыдущие исследования касались фазовой диаграммы, РТ-условий диссоциации и определения стабильных фаз бензола при давлениях до 30 ГПа (Thiery, Leger, 1988; Ciabini et al., 2005; 2007; Kondrin et al., 2017; Chanyshev et al., 2018). Теоретический расчет фазовых преобразований бензола проведен до давлений 200 ГПа (Wen et al., 2011). В результате моделирования были установлены параметры полимеризации бензола в графеноподобную структуру, а также узкий интервал стабильности проводящей фазы в районе 190 ГПа. Экспериментальное исследование P-Vуравнений состояния проводились для нафталина (Fabbiani et al., 2006; Likhacheva et al., 2014a; 2014b), антрацена, тетрацена, пентацена (Oehzelt et al., 2006) и коронена (Zhao et al., 2013; Чанышев и др., 2016). Теоретические расчеты уравнений состояния первопринципными методами проводились в отдельных работах практически для всех рассмотренных ПАУ, за исключением фенантрена, аценафтена и пирена. Результаты основные работ (Venuti et al., 2002; 2004; Oehzelt et al., 2006; Журавлев и др., 2012)обсуждаются ниже. Также для некоторых фаз экспериментально установлены высокобарические модификации, однако корректное определение и описание структур для них пока не проводилось (Zhao et al., 2013).

В данной работе впервые систематизированы параметры сжимаемости ПАУ в зависимости от их молекулярной и кристаллической структуры, а также от размера молекулы и количества атомов в элементарной ячейке. Это имеет критическое значение для их дальнейшего экспериментального и теоретического исследования при высоких давлениях и температурах.

МЕТОДИКА РАСЧЕТОВ

Параметры сжимаемости ПАУ были получены на основании расчетов электронной структуры, выполненных в рамках теории функционала плотности, который является одним из наиболее широко используемых методов вычислительной физики и квантовой химии (Dreizler, Gross, 1995). Основным преимуществом данного метода является то, что он позволяет существенно упростить задачу описания многоэлектронной системы путем перехода от рассмотрения многоэлектронной волновой функции, зависящей от 3N_e переменных (N_e – число электронов в системе), к исследованию функции распределения электронной плотности, которая в пренебрежении электронным спином зависит только от трех пространственных координат (Hohenberg, Kohn, 1964).

Расчеты методом функционала плотности проводятся в рамках формализма Кона-Шэма (Kohn, Sham, 1965), в котором сложная задача описания системы сильновзаимодействующих электронов и ядер сводится к более простой задаче. В ней рассматривается система невзаимодействующих электронов, движущихся в эффективном потенциале, в котором помимо электростатического потенциала ядер учитывается кулоновская, обменная и корреляционная энергия электронов. Описание последних двух взаимодействий представляет основную трудность метода теории функционала плотности, и различные попытки их учесть объясняют большое разнообразие функционалов, разработанных к настоящему времени. Другой трудностью метода является то, что дисперсионные взаимодействия не описываются им с должной точностью. Для преодоления этой трудности разработаны полуэмпирические поправки к теории, которые применяются и в настоящей работе.

Расчеты выполнены в рамках теории функционала плотности в программном пакете VASP (Kresse, Furthmüller, 1996; Kresse, Joubert, 1999) с градиентной аппроксимацией обменно-корреляционного потенциала в форме optPBE-vdW (Klimeš et al., 2010). Для корректного учета ван-дер-Ваальсовых взаимодействий, вклад которых существенен в исследуемых молекулярных соединениях, использовался метод vdW-DF (Dion et al., 2004; Klimeš et al., 2011), реализованный в программе VASP по схеме, предложенной в работе (Román-Pérez, Soler, 2009). Данный метод был нами успешно применен ранее при исследовании свойств кристалла нафталина (Likhacheva et al., 2014b). Энергия обрезания плоских волн была равна 600 эВ. Для расчета электронной структуры использовали необходимое (в зависимости от вещества) число k-точек в зоне Бриллюэна по схеме Монкхорста-Пака, равное примерно 2000/N, где N – число атомов в элементарной ячейке. Расчет проводился для одной элементарной ячейки каждого вещества с 2 или 4 молекулами на ячейку в зависимости от структуры. Координаты атомов в элементарной ячейке были оптимизированы с использованием метода сопряженных градиентов. Форма элементарной ячейки оптимизировалась при условии сохранения ее объема.

Для определения границ абсолютной устойчивости кристаллов использовалось условие ∂P/∂V = 0. Устойчивость твердых тел должна сопровождаться выполнением очевидного условия ∂P/∂V < 0, означающего то, что при сжатии давление в системе растет. Нарушение данного условия приводит к неустойчивости кристаллической структуры и свидетельствует о высокой вероятности фазового перехода, полимеризации или аморфизации исходного вещества.

Для расчета упругих параметров ПАУ полученные данные интегрировали по уравнению P = −∂F/∂V и далее аппроксимировали, используя уравнение состояния Винэ (Vinet et al., 1987)в интервале давлений 0–20 ГПа:

$P = 3K_0 \frac{1 – x^{{\displaystyle \frac{1}{3}}}}{x^{{\displaystyle \frac{2}{3}}}} exp\left[\eta \left(1–x^{\frac{1}{3}}\right)\right],$

где $х=V/V_0$ и $\eta =\frac{3}{2(K\text{'}–1)}$, K – модуль всестороннего сжатия и K´ – его производная по давлению. Методика расчета идентична той, которую мы применяли для исследования нафталина (Likhacheva et al., 2014b).

ПОЛИМОРФНЫЕ ПЕРЕХОДЫ В БЕНЗОЛЕ И ПАУ ПРИ ВЫСОКОМ ДАВЛЕНИИ

Набор исследованных соединений представлен на рис. 1 и в табл. 1. Были выбраны ПАУ от нафталина до коронена с учетом структуры молекулы для соединений из трех, четырех и шести бензольных колец. Теоретический расчет стабильных структур ПАУ при высоких давлениях является отдельной задачей. Пока сложно проводить свободный поиск структур для веществ с большим количеством атомов в элементарной ячейке и, кроме того, среди молекулярных соединений, где существенным является подбор псевдопотенциалов, учитывающих Ван-дер-Ваальсовы взаимодействия. Поэтому в данной работе мы ориентировались на структуры, стабильные при нормальных условиях, – 298.15 K и 0.1 МПа. Для большинства соединений характерна структура P2₁ /с или Р2₁ /а или сходная с ней P-1. Эти структуры весьма стабильны при повышении давления. Например, подобная структура бензола термодинамически стабильна до давления 190 ГПа (Wen et al., 2011). Тем не менее для ряда соединений установлены высокобарические модификации, суммированные в табл. 1. Следует отметить, что часть фазовых переходов, установленных для ПАУ при высоких давлениях, в последствии не подтверждается, а часть полиморфных переходов происходит внутри одной пространственной группы (Venuti et al., 2004; Fabbiani et al., 2006). Весьма вероятным видится появление метастабильных или полимерных фаз за счет сдвиговых напряжений в эксперименте (Политов и др., 2000). Полимеризация или аморфизация ПАУ наблюдается во многих экспериментах, при этом давление подобных переходов может сильно варьировать (Zhao et al., 2013; Чанышев и др., 2016).

Рис. 1. Молекулярная структура исследованных бензола и ПАУ.

Таблица 1. Известные структуры и параметры фазовых переходов бензола и ПАУ при 298 K и высоком давлении

Бензол имеет сложную фазовую диаграмму, исследованную экспериментально до давлений около 30 ГПа (Thiery, Leger, 1988; Ciabini et al., 2005; 2007). Теоретический поиск стабильных структур проведен при давлении до 200 ГПа (Wen et al., 2011). Экспериментальные и теоретические данные хорошо согласуются в области умеренных давлений. Номенклатура фаз варьирует в разных работах, мы придерживались работы (Wen et al., 2011), где приняты следующие обозначения: I – Pbca, II – P4₃2₁2, III – P2₁ /c (P2₁/a) (табл. 1). Согласно экспериментальным данным, фаза III аморфизуется при давлении около 40 ГПа на 298 K. Теоретические расчеты показывают стабильность фазы III до 40–50  ГПа, где она превращается в фазу V (P2₁) (Wen et al., 2011). Высокотемпературная фаза IV(стабильная при 3–6 ГПа и 600–800 К) была предложена в работе (Cansell et al., 1993). Однако более детальные исследования фазовой диаграммы бензола с помощью многопуансонной техники не подтвердили наличие этой фазы (Chanyshev et al., 2018). Следует отметить возможную полимеризацию бензола при высоком давлении и образование графеноподобных структур (Wen et al., 2011), которые могут быть характерны и для всех ПАУ.

Для пирена и фенантрена определены высокобарические модификации со структурами, близкими к стабильным при 0.1 МПа. При этом сильно меняются параметры ячейки без существенного изменения кристаллической структуры (Fabbiani et al., 2006). Для нафталина, антрацена, тетрацена и пентацена исследована сжимаемость до давлений 9–12 ГПа (23 ГПа для антрацена) – при этом не обнаружено фазовых переходов (табл. 1). Однако более детальный анализ экспериментальных данных, а также последующие теоретические расчеты показали полиморфные переходы в тетрацене и пентацене при 5–7 ГПа без изменения структуры (Schatschneider et al., 2013). Фазовые соотношения при высоких давлениях не определялись для аценафтена и гексацена.

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

Результаты первопринципных расчетов представлены в табл. 2 и 3. В табл. 2 проводится сравнение параметров элементарной ячейки при нулевом давлении с предыдущими экспериментальными и теоретическими данными. Учитывая примерное значение теплового расширения для молекулярных соединений (α~10^-4  К^-1), различия между рассчитанными и измеренными параметрами находятся в пределах 3%, что свидетельствует о надежности проведенных расчетов. Сравнение сжимаемости с экспериментальными и предыдущими теоретическими расчетами проведено на рис. 2, 3. Оно также показывает близкое соответствие между расчетами и экспериментами.

Таблица 2. Рассчитанные параметры элементарной ячейки при стандартных условиях в сравнении с экспериментальными данными и предыдущими теоретическими расчетами

В ходе расчетов для всех рассмотренных кристаллических решеток наблюдается их хорошая гидростатичность, которая улучшается с увеличением давления. Эти данные были получены в результате анализа рассчитанных величин тензора напряжения, диагональные члены которого в случае идеального гидростатического сжатия должны быть равны между собой, а недиагональные – обращаться в ноль. Так, абсолютные значения недиагональных компонент тензора напряжений не превышают величин 0.05 ГПа для всех изученных соединений во всем диапазоне давлений. Разница между ними достигает величин 0.05–0.07 ГПа, что составляет в некоторых случаях 70% от абсолютных величин компонент тензора в интервале давлений ниже 0.1 ГПа. Однако с ростом давления относительные величины отклонения диагональных членов тензора напряжений резко уменьшаются и не превышают 5% при давлениях до 3 ГПа и 1% при более высоких давлениях.

Рис. 2. Сравнение сжимаемости бензола III с литературными данными. 1 – данная работа, 2 – аппроксимация по УС Винэ, 3 – экспериментальные данные (Ciabini et al., 2005) при 540 К и их аппроксимация по УС Винэ (4), 5 – экспериментальные данные (Thiery, Leger, 1988) при 294 К, 6 – расчет из первых принципов (Wen et al., 2011).

Рис. 3. Сравнение сжимаемости ПАУ, полученной в данной работе, с литературными данными. (а) Нафталин: 1 – данная работа и аппроксимация по УС Винэ (2), 3 – экспериментальные данные (Likhacheva et al., 2014a; 2014b) и их аппроксимация по УС Винэ (4), 5 – расчет из первых принципов (Журавлев и др., 2012). (б) Антрацен: 1 – данная работа и аппроксимация по УС Винэ (2), 3 – экспериментальные данные и 4 – расчет из первых принципов (Oehzelt et al., 2006), 5 – расчет из первых принципов (Журавлев et al., 2012). (в) Пентацен: 1 – данная работа и аппроксимация по УС Винэ (2), 3 – экспериментальные данные и 4 – расчет из первых принципов (Oehzelt et al., 2006).

На рис. 2 показаны данные для бензола III. Следует отметить плохую сходимость между экспериментальными измерениями в работах (Thiery, Leger, 1988) при 295 K и (Ciabini et al., 2005) при 540 K. Данные при 540 K соответствуют меньшим объемам, чем данные при 295 K. Это может объясняться тем, что температура снижает влияние девиаторных напряжений, которые влияют на точное определение объема ячейки. Кристаллы в образце в алмазной ячейке могут ориентироваться вдоль самого жесткого направления параллельно плоскостям алмаза, и это может повлиять на измерения. Очевидно, что нужны дополнительные исследования сжимаемости бензола. Наши расчетные данные в целом согласуются с объемами, измеренными в работе (Ciabini et al., 2005), однако параметры сжимаемости (табл. 3) также близко соответствуют экспериментальным данным из работы (Thiery, Leger, 1988).

Таблица 3. Значения параметров упругости по УС Винэ в интервале давлений 0–20 ГПа в сравнении с предыдущими экспериментальными и расчетными данными

Примечания. A – данная работа; эксп. – экспериментальные данные; теор. – теоретический расчет. 1 – (Ciabini et al., 2005); 2 – (Likhacheva et al., 2014b); 3 – (Журавлев и др., 2012); 4 – (Oehzelt et al., 2006). В скобках показано стандартное отклонение для последних значащих цифр.

Рис. 4. Аномальная зависимость параметра а пентацена от давления. Черные кружки – данная работа, белые – данные работы (Schatschneider et al., 2013).

Сравнение с экспериментальными данными для нафталина, тетрацена и пентацена также показывает близкое соответствие полученных упругих параметров экспериментальным значениям (рис. 3). Максимальные различия для пентацена составляют 4% при 20 ГПа. Однако пентацен, как и аценафтен, имеет аномальную сжимаемость с возможными фазовыми переходами, которые следуют из расчетов. Аценафтен со структурой Pmc2₁ имеет границу абсолютной устойчивости при давлении около 1.5 ГПа, поэтому расчет для него показывает условные параметры сжимаемости.

Пентацен имеет аномалии зависимости параметров элементарной ячейки от давления при 6–7 ГПа (рис. 4) и соответствующее небольшое изменение объема ячейки (рис. 3). Это уже было показано в работе (Schatschneider et al., 2013), где отмечается, что подобная зависимость связана с полиморфным переходом. При этом не меняются симметрия и мотив кристаллической решетки, а происходит лишь увеличение плотности упаковки за счет увеличения межплоскостных углов или угла между молекулами (елочного угла, herringbone angle). Там же показано, что наиболее заметны эти изменения при расчете поверхности Хиршфельда (Прокаева и др., 2009; Schatschneider et al., 2011). В работе (Schatschneider et al., 2013) подобные аномалии установлены и для тетрацена при давлении около 5 ГПа, однако в наших расчетах они не были обнаружены.

Расчет для различных полиморфов сделан для бензола и фенантрена. Он показывает, что сжимаемость двух различных полиморфов практически не отличается (табл. 3, рис. 5). Поэтому сравнение сжимаемостей между различными ПАУ правомерно, с нашей точки зрения, даже для структур, стабильных при нулевом (атмосферном) давлении, если структура не теряет абсолютной стабильности, которая фиксируется по выполнению условия ∂P/∂V = 0 при высоком давлении. Безусловно, последующие расчеты, теоретический поиск стабильных структур или новые экспериментальные определения полиморфных переходов дадут более точные результаты подобного сравнения.

Рис. 5. Сравнение кривых сжимаемости бензола и ПАУ, полученных в данной работе: (а) на интервале давлений 0–25 ГПа, (б) в узком интервале давлений, позволяющем увидеть различия между отдельными линиями.

Рис. 6. Зависимость модуля всестороннего сжатия (K₀) от молекулярной массы ПАУ. Пунктирной линией показан тренд увеличения K₀.

Все кривые сжимаемости исследованных ПАУ показаны на рис. 5. Очевидно, что все вещества имеют довольно близкие параметры сжимаемости (K₀ и K’) с небольшим отклонением для бензола (табл. 3, рис. 5). Разница между относительным объемом всех ПАУ от нафталина до коронена при 20 ГПа не превышает 5%. Производные K’ находятся в узком интервале – 6.9–7.5, поэтому можно проводить корреляцию модулей сжатия без дополнительной нормировки на единое значение K’. Наблюдается выраженная зависимость в ряду от бензола до гексацена (с увеличением атомной массы вещества) в сторону увеличения K₀ (рис. 6). Из общей зависимости выбиваются лишь аценафтен и пентацен. Скорее всего, это связано с отмеченными выше аномалиями сжимаемости из-за полиморфных переходов. Не наблюдается зависимости сжимаемости от кристаллической структуры. Как отмечалось, все вещества имеют близкую структуру, несмотря на различия в пространственных группах. Среди веществ с одинаковым числом бензольных колец большую сжимаемость имеют менее плотные структуры (например, тетрацен и гексацен по сравнению с пиреном и короненом, рис. 5).

Полученные данные по сжимаемости ПАУ можно использовать при разработке высокотемпературных уравнений состояния и расчета равновесного состава жидких и твердых компонентов системы С-О-Н. Сложность моделирования ПАУ при высоких температурах заключается в их последовательной полимеризации с увеличением температуры (Чанышев и др., 2014a;
Chanyshev et al., 2015b; 2017; Kondrin et al., 2017). Это обуславливает различия в полях стабильности, установленных экспериментально (Davydov et al., 2004; Чанышев и др., 2014a; Chanyshev et al., 2015a) и в результате термодинамического моделирования (Карпов и др., 1998; Зубков, 2001; 2009).

ЗАКЛЮЧЕНИЕ

В работе приведены результаты первопринципных расчетов кристаллической структуры и параметров сжимаемости бензола и ряда полициклических ароматических углеводородов. Использовался метод теории функционала плотности с градиентной аппроксимацией обменно-корреляционного потенциала в форме PBEс учетом Ван-дер-Ваальсовых взаимодействий (optPBE-vdW) с помощью программы VASP в интервале давлений 0–20 ГПа. Сравнение с экспериментальными данными и предыдущими теоретическими расчетами показывает релевантность выбранных подходов.

Все исследованные вещества от бензола до коронена и гексацена (с шестью бензольными кольцами) имеют близкую сжимаемость с модулем всестороннего сжатия от 8 до 12 ГПа и его производной по давлению 6.9–7.5. Близость рассчитанных параметров указывает на основную роль уменьшения межмолекулярных расстояний при сжатии ПАУ и слабую деформацию самих молекул и бензольных колец. Установлена слабая зависимость сжимаемости от количества атомов (бензольных колец) в молекуле вещества или от типа кристаллической структуры. Соединения с большим количеством бензольных колец, а также более плотной структурой колец – циклические (пирен, коронен), имеют меньшую сжимаемость по отношению к менее плотным ПАУ (тетрацен, гексацен). Некоторые ПАУ имеют высокобарические модификации, однако корректное описание их структур, позволяющее получить упругие модули, еще не сделано и требует дополнительных исследований.

Благодарности

Авторы выражают благодарность Информационно-вычислительному центру Новосибирского государственного университета за возможность использования вычислительного кластера, на котором были произведены все расчеты.

Источник финансирования

Работа выполнена по государственному заданию ИГМ СО РАН.

About the authors

K. D. Litasov

Author for correspondence.
Email: klitasov@igm.nsc.ru
Russian Federation, 3 Koptyuga Ave., Novosibirsk; 2 Pirogova St., Novosibirsk

T. M. Inerbaev

Email: talgat.inerbaev@gmail.com
Russian Federation, 3 Koptyuga Ave., Novosibirsk; 2, Satpayev str., Astana

F. U. Abuova

Email: super.fatima-82@yandex.kz
Kazakhstan, 2, Satpayev str., Astana

A. D. Chanyshev

Email: chanyshev_90@mail.ru
Russian Federation, 3 Koptyuga Ave., Novosibirsk; 2 Pirogova St., Novosibirsk

A. K. Dauletbekova

Email: dauletbekova_ak@enu.kz
Kazakhstan, 2, Satpayev str., Astana

A. T. Akilbekov

Email: akilbekov_at@enu.kz
Kazakhstan, 2, Satpayev str., Astana

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Molecular structure of the studied benzene and PAH.

Download (160KB)

Indexing metadata

3. Fig. 2. Comparison of the compressibility of benzene III with the literature data. 1 - this work, 2 - approximation by the Vine CS, 3 - experimental data (Ciabini et al., 2005) at 540 K and their approximation by the Vine CSS (4), 5 - experimental data (Thiery, Leger, 1988) at 294 K, 6 is a first-principle calculation (Wen et al., 2011).

Download (104KB)

Indexing metadata

4. Fig. 3. Comparison of the compressibility of PAHs obtained in this work with the literature data. (a) Naphthalene: 1 - this work and approximation by the CSS Vine (2), 3 - experimental data (Likhacheva et al., 2014a; 2014b) and their approximation by the CSS Vine (4), 5 - calculation from first principles (Zhuravlev et al., 2012). (b) Anthracene: 1 - this work and approximation by the CSS Vine (2), 3 - experimental data and 4 - calculation from first principles (Oehzelt et al., 2006), 5 - calculation from first principles (Zhuravlev et al., 2012). (c) Pentacene: 1 — this work and approximation by the CSS Vine (2), 3 — experimental data and 4 — first-principle calculations (Oehzelt et al., 2006).  

Download (282KB)

Indexing metadata

5. Fig. 4. Anomalous dependence of the parameter a pentacene on pressure. The black circles are the given work, the white ones are the data of the work (Schatschneider et al., 2013).  

Download (96KB)

Indexing metadata

6. Fig. 5. Comparison of compressibility curves of benzene and PAHs obtained in this work: (a) at a pressure range of 0–25 GPa, (b) in a narrow pressure range, which allows to see the differences between individual lines.  

Download (353KB)

Indexing metadata

7. Fig. 6. Dependence of the bulk compression modulus (K0) on the molecular weight of PAHs. The dotted line shows the trend of increasing K0.

Download (88KB)

Indexing metadata