Журнал физической химии, 2020, T. 94, № 6, стр. 855-861

ДЕТАЛИ РАСЧЕТОВ

Методом молекулярной динамики в режиме NPT были получены траектории, усредненные за 0.1 пс для получения колебательно-усредненных V-структур 1,3-ПД при трех температурах 293, 323 и 373 K и давлении 1 бар. Модель жидкого 1,3-ПД с 1000 молекулами в прямоугольной ячейке с периодическими граничными условиями была построена методом молекулярной динамики [8, 9]. Параметры взаимодействия и парциальные заряды на атомах были такие же, как и для ЭГ [8] (подробности см. в [9]). Для численного интегрирования уравнений движения использовался скоростной алгоритм Верле [10]. Шаг интегрирования составлял 1 фс. Температура в системе поддерживалась с помощью столкновительного термоcтата [11, 12], давление – с помощью баростата Берендсена [13].

При каждом значении температуры для изучаемых систем проводились предварительные расчеты их молекулярной динамики в режиме NPT до достижения стационарного поведения, при котором характеристики систем выходили на плато и флуктуировали около своих средних значений. Длина таких траекторий достигала 1 нс. Полученные при температурах 293, 323 и 373 K средние плотности ЭГ были равны 1.083, 1.056 и 1.016 г/см³ соответственно; плотности 1,3-ПД – 1.088, 1.066 и 1.024 г/см³ соответственно.

Из полученных методом молекулярной динамики функций радиального распределения был определен геометрический критерий водородной связи (как межмолекулярной, так и внутримолекулярной). Для конструирования такого критерия были использованы положения первых минимумов этих функций для атомов кислорода и водорода гидроксильных групп. Наличие водородной связи регистрировалось, если расстояние между кислородами двух гидроксильных групп не превышало 3.5 Å, а расстояние между кислородом одной из гидроксильных групп и протоном другой гидроксильной группы не превышало 2.5 Å.

Усреднение траекторий жидкой системы по небольшому интервалу времени сглаживает движение атомов, устраняя их тепловые высокочастотные флуктуационные колебания около текущих средних положений и, тем самым, исключает регистрацию согласно геометрическому критерию “нефизических” водородных связей, вызванных такими флуктуациями. V-структуры ЭГ и выбор временного интервала усреднения обсуждались в работе [14]. Взятый здесь интервал усреднения 0.1 пс устраняет высокочастотные флуктуационные движения и валентные колебания атомов и одновременно учитывает наличие более медленных структурных изменений в жидкости, приводящих к образованию и разрыву водородных связей.

Базируясь на критерии Н-связи, построили матрицу смежности для V-траекторий, полученных усреднением координат атомов за время 0.1 пс. В каждый момент времени V-траектория содержит структуру из 1000 молекул в расчетной ячейке с периодическими граничными условиями. Элемент с индексом (i, j) матрицы смежности содержит 1, если i-я молекула связана с j-й молекулой водородной связью (удовлетворяет критерию Н-связи), и 0, если молекулы не связаны; то же относится и к внутримолекулярной Н-связи.

Основным параметром, который оценивает сетку водородных связей внутри структуры, служит размер максимального компонента связности графа, описываемого построенной матрицей смежности. Для нахождения компонент связности графа был использован алгоритм Тарьяна [15], который анализирует матрицу смежности. Результатом его работы является количество компонент связности в графе и их размер, а также список вершин графа с указанием их принадлежности к тому или иному компоненту связности.

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

Результаты численного анализа V-траекторий при температурах 293 и 373 K для 1,3-ПД и этиленгликоля приведены в табл. 2.

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

Основной результат анализа табл. 2 – вывод о том, что практически все молекулы во всех исследованных траекториях 1,3-ПД и ЭГ образуют трехмерную сетку водородных связей, только при повышенной температуре минимальный размер максимального связного компонента у ЭГ равен 996, а у 1,3-ПД – 992, что соответствует отрыву четырех и восьми молекул. Отметим большее количество внутримолекулярных Н-связей у 1,3-ПД, увеличивающееся с возрастанием температуры, что согласуется с данными ИК-спектроскопии [2].

Важнейшая характеристика пространственной сетки Н-связей – время жизни водородной связи. Рассмотрим суммарное время жизни водородных связей в 1,3-ПД и ЭГ. Для этого построим матрицы смежности в каждый временной отсчет по 0.1 пс и просуммируем их вдоль всей траектории V-структур. Основные характеристики полученной матрицы для обоих веществ при 293, 323 и 373 K приведены в табл. 3; приводятся также характеристики межмолекулярных связей для структур 1,3-ПД и ЭГ, в каждой структуре 1000 молекул.

Количество пар молекул, образовывавших водородную связь, – это суммарное количество пар молекул, между которыми хотя бы один раз вдоль траектории оказался выполненным критерий водородной связи. Число Н-связей, приходящихся на одну молекулу, усредненное по структуре. Как и ранее в табл. 2, приведены три значения: первое – усредненное значение по всем структурам, а в квадратных скобках приведены минимальное и максимальное значение параметра по всем структурам траектории.

Время жизни водородной связи вычислялось следующим образом. Если критерий водородной связи оказался выполненным для некоторой пары молекул в структуре, то, поскольку V-траектории усреднены за время 0.1 пс, мы считаем, что критерий был выполнен в течение 0.1 пс. Просуммировав для выбранной пары молекул время выполнения критерия водородной связи вдоль всей траектории, получаем интегральное время жизни данной водородной связи. На рис. 1 для 1,3-ПД и ЭГ при температурах 273 и 373 K приведены гистограммы распределения долей числа пар молекул, образовывавших на траектории Н-связи с интегральными временами их жизни, попадающими в заданные интервалы времени. Количество молекул, взаимодействующих с данной молекулой на траектории заданной длины, – это суммарное число других молекул жидкости, хотя бы один раз образовавших водородную связь с выбранной молекулой.

Рис. 1.

Гистограммы распределения долей числа пар молекул, образовывавших на траектории Н-связи с интегральным временем их жизни попадающим в соответствующий временной интервал для 1,3-ПД (а) и ЭГ (б) при 373 (закрашенные области) и 273 K (контуры). Длина интервала – 4 пс.

Длительности всех траекторий немного превышают 1 нс, и можно считать, что они различаются незначительно, так что рассматриваемые суммарные показатели сравнимы между собой. При температуре 293 K у ЭГ количество пар молекул, образовывавших водородную связь в ∼1.55 раза больше, чем у 1,3-ПД. С ростом температуры соотношение меняется, и при 373 K у ЭГ таких пар молекул становится в _˜1.9 раза больше, чем у 1,3-ПД. Число водородных связей, приходящихся на одну молекулу с ростом температуры у обоих веществ немного уменьшается, при этом само значение показателя немного больше у ЭГ, чем у 1,3-ПД, и с ростом температуры оно не меняется.

Время жизни водородной связи при повышении температуры сильно падает. Для ЭГ этот показатель уменьшается более, чем в 5 раз, а для 1,3‑ПД – в 4 раза. Максимальное время жизни Н-связи при 293 K, например, у 1,3-ПД составляет 1088.3 пс. Это означает, что имеются пары молекул с неразрывной водородной связью практически вдоль всей траектории. Средние значения времен жизни водородной связи у 1,3-ПД в ∼1.5 раза превышают соответствующие значения у ЭГ. Интегральное время жизни водородной связи характеризует время, через которое две близко расположенные молекулы в значительной степени меняют свое окружение, так, что их новые контакты становятся маловероятными. Кажется, что времена жизни межмолекулярных водородных связей должны коррелировать с вязкостью. Такое соответствие действительно есть, но только качественное. Так, с увеличением температуры от 293 до 373 K среднее интегральное время жизни водородной связи для ЭГ монотонно уменьшается немного более, чем в 5 раз, в то время как экспериментально определенная вязкость уменьшается в _˜10 раз. Еще меньшая корреляция наблюдается с непрерывным временем жизни Н-связи. По-видимому, связь между временами жизни водородных связей и вязкостью имеет более сложную природу и требует проведения дополнительных исследований.

Для определения времени жизни конкретной водородной связи был выбран подход Волошина [16]. В соответствии с этим подходом и выбранным критерием водородной связи были получены характеристики непрерывного времени жизни Н-связи, времени до первого разрыва относительно выбранного нами момента времени t_о, как и в продолжении его (t_o + Δt), так и в обратном направлении (t_o – Δt). В наших расчетах Δt = 200 пс. Важно подчеркнуть, что речь идет об анализе всех водородных связей, зарегистрированных в заданный момент времени (t_o). Непрерывное время жизни каждой Н-связи определяется суммой времен ее жизни без прерываний до этого момента и после него. Это непрерывное время жизни Н‑связи положим в основу перколяционного анализа пространственных сеток для выделенных конфигураций. Под конфигурацией понимается структура (набор водородных связей) в момент времени t_o со специальной матрицей смежности. Методика расчета представлена в [7]. Будем называть ребра окрашенными, если время жизни этой Н-связи не меньше некоторого времени t*. Рассмотрим подграф, в который входят только окрашенные связи, на условие “протекания”, т.е. образует ли этот подграф трехмерную сетку Н-связей с учетом периодических граничных условий для нашей структуры. Если условие “протекания” выполняется, то самое большое из возможных значений $t_{{{\text{кр}}}}^{*}$ будем называть порогом перколяции. Связный подграф, образующий трехмерную сетку на пороге перколяции, называется перколяционным кластером.

Результаты расчетов порога перколяции и мощности перколяционного кластера при 293, 323 и 373 K для 1,3-ПД и ЭГ приведены в табл. 4. Перколяционный анализ проводился с усреднением результатов по 10 независимо выбранным моментам времени t_o внутри траектории и конфигурациям, соответствующим этим временам.

В табл. 4 приведены следующие характеристики рассмотренных конфигураций:

– среднее число связей на молекулу (рассчитывается как отношение общего количества Н-связей к числу молекул в конфигурации); максимальное время жизни непрерывной связи в конфигурации – самое большое значение времени жизни непрерывной связи, полученное в данной конфигурации;

– среднее время жизни непрерывной связи (усредненное по всем водородным связям в конфигурации значение);

– порог перколяции (время жизни водородной связи, при котором подграф рассматриваемой конфигурации, содержащий водородные связи с временем непрерывной жизни не меньше, чем время образования трехмерной бесконечной сетки);

– среднее время жизни окрашенных связей на пороге протекания (усредненное значение времени непрерывной жизни водородных связей, если рассматривать только те связи, которые существовали по времени дольше, чем порог протекания);

– доля окрашенных связей на пороге протекания (отношение количества связей, которые существовали дольше, чем порог протекания, к общему количеству связей в конфигурации);

– доля молекул в перколяционном кластере (отношение количества молекул, входящих в перколяционный кластер, к общему количеству молекул);

– доля окрашенных связей в перколяционном кластере (показатель, аналогичный доле окрашенных связей на пороге протекания, только при его расчете рассматриваются лишь те водородные связи, которые входят в перколяционный кластер; то же самое относится и к паре показателей: среднее время жизни окрашенных связей в перколяционном кластере и среднее время жизни окрашенных связей на пороге протекания.

Среднее число связей на молекулу для обоих веществ уменьшилось на _˜ 10% с ростом температуры. Этот показатель для 1,3-ПД на те же 10% меньше, чем для этиленгликоля при всех рассмотренных температурах. Все показатели времени жизни в табл. 4 с ростом температуры сильно уменьшаются: и максимальное время жизни непрерывной связи в конфигурации, и среднее время жизни непрерывной связи, и среднее время жизни окрашенных связей на пороге протекания, и среднее время жизни окрашенных связей в перколяционном кластере. Однако, по сравнению с уменьшением времени жизни водородной связи вдоль всей траектории, приведенным в табл. 3, скорость падения у этиленгликоля несколько меньше, чем у 1,3-ПД. Порог перколяции этиленгликоля при 373 K в 2.9 раза меньше, чем при 293 K. Такой же параметр для 1,3-ПД уменьшился в 3.625 раза. Все временные характеристики для 1,3-ПД несколько больше, чем для этиленгликоля.

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

Таким образом, описание колебательно-усредненных V-структур ЭГ и 1,3-ПД, полученных методом молекулярной динамики при трех температурах 293, 323 и 373 K, проведенное методами теории графов и перколяции, показало наличие трехмерной сетки водородных связей, в которую включены почти все молекулы рассматриваемых структур. Получены характеристики структур. Отмечено, что важнейшей характеристикой структуры является время жизни водородной связи. По всем временным характеристикам водородная связь в 1,3-ПД более долгоживущая, чем в структуре ЭГ, что объясняет более широкий температурный интервал жидкой фазы и большую вязкость 1,3-ПД по сравнению с ЭГ. Отмечено большее количество внутримолекулярных Н-связей у 1,3-ПД, увеличивающееся с возрастанием температуры. Несмотря на большее количество –СН₂-групп в молекуле 1,3-ПД, которые рассматриваются как дефекты сетки [6], обе сетки (у ЭГ и 1,3-ПД) устойчивы и похожи, возможно, потому, что молекула 1,3-ПД более гибкая, чем молекула ЭГ, и именно из-за того, что в ее молекуле не две, как у ЭГ, а три –СН₂-группы.

Применение топологических методов вместе с методом молекулярной динамики к описанию пространственных сеток Н-связей в жидкостях и растворах позволяет получать пространственные и временные характеристики этих сеток, объясняющие физико-химические свойства жидкостей, обладающих такими сетками, характеристики, которые нельзя получить прямым экспериментом.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта № 19-03-00215). Основные расчеты проведены в МСЦ РАН и ИПМ имени М.В. Келдыша РАН.