Журнал физической химии, 2022, T. 96, № 1, стр. 86-95

Моделирование растворов перхлората лития в сульфолане методом классической молекулярной динамики

Моделирование растворов перхлората лития в сульфолане методом классической молекулярной динамики проводили в программном пакете Gromacs version 2019.1 [16]; использовали силовые поля, входящие в программный пакет GROMACS: Amber (GAFF) [17, 18], Oplsaa [19], Gromos [20] и Charmm [21].

Недостающие параметры (длины связей, валентные и тетраэдрические углы, заряды) для молекул сульфолана и аниона перхлората рассчитывали в программном пакете Gaussian 09 [22] методом B3LYP/ aug-cc-pVTZ. Выходные файлы с оптимизированным строением молекул или ионов преобразовывали в файлы, поддерживаемые программой Gromacs, с помощью интернет сервисов [23–29]. Парциальные заряды на атомных центрах молекул и анионов рассчитывали методами Маллекена и Мерц–Коллмана. Также парциальные заряды атомных единиц молекул сульфолана брали из [23, 30, 31].

Константы, описывающие силовые поля, рассчитывали с помощью интернет-сервисов [23–29]. Параметры силовых полей приведены во вспомогательном материале.

Подготовку системы к МД проводили согласно процедуре, описанной в [32]. В бокс размером 50 Å × 50 Å × 50 Å с периодическими граничными условиями помещали заданное количество ионов фоновой соли и/или молекул растворителя. Соотношение ионов фоновой соли и молекул растворителя было 1/10 для моделирования 1 М раствора LiClO₄ в сульфолане. После минимизации энергии с критериями сходимости 1000 кДж/моль проводили уравновешивание системы следующим образом. Первоначально уравновешивали систему в изобарно-изотермическом ансамбле (NPT) в течение 400 пс при 200 атм и 303 К, затем 20 нс при 1 атм и 303 К.

Для расчета плотности системы выполняли симуляцию в NPT-ансамбле для сбора данных 20 нс при 1 атм и 303 К. Для расчета вязкости, коэффициентов самодиффузии, электропроводности выполняли симуляцию в изохорно-изотермическом (NVT) ансамбле в течение 10–40 нс.

Время интегрирования составляло 1 фс. Постоянную температуру (Т = 303 К) поддерживали термостатом Нозе–Хувера [33, 34] с параметром τ = 0.2 пс. Давление поддерживали постоянным (1 атм) баростатом Берендсен [35] с параметром τ = 0.2 пс.

В начальный момент времени атомам придавали случайные (по величине и направлению) скорости, нормированные таким образом, чтобы кинетическая энергия системы соответствовала заданной температуре. Для каждого объекта выполняли несколько последовательных моделирований для подсчета данных и оценки погрешностей. В последовательных моделированиях результаты – статистически независимы и нормально распределены. Погрешность расчета плотности не превышала 0.5%, а электропроводности, коэффициентов диффузии и вязкости – 15%.

Размеры молекул сульфола и аниона перхлората лития рассчитывали в программном пакете Chemcraft version 1.8. Диаметр частиц рассчитывали как диаметр столкновения частиц (collision diameter). Объекты моделирования представлены на рис. 1.

Рис. 1.

Объекты моделирования: молекула сульфолана (а), анион перхлората (б) и катион лития (в).

Анализ данных. Расчет физико-химических свойств растворов

Коэффициенты диффузии рассчитывали с использованием соотношения Эйнштейна [36] из среднеквадратичного смещения атомов из набора начальных положений с помощью пакетов моделирования, реализованных в GROMACS.

Электропроводность рассчитывали двумя методами – из автокорреляционной функции тока (соотношение вращательного и поступательного дипольного момента системы) с помощью аппроксимации Эйнштейна–Гельфанда и по уравнению Нернст–Эйнштейна [3, 37].

Вязкость рассчитывали из корреляционной функции поперечного тока с использованием полей поперечного импульса (автокорреляционная функция поперечного тока) с помощью пакетов моделирования, реализованных в GROMACS. Всего учитывали 16 автокорреляционных функций поперечного тока, соответствующих различным k-векторам, что дает 16 значений η. Значения η соответствуют η(k) = η_∞(1 – ak²), что дает сдвиговую вязкость η_∞ при k = 0.

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

Координационное число молекул типа i в первичной сольватной оболочке, окружающей одиночную частицу j, рассчитывали интегрированием радиальной функции распределения с помощью пакетов моделирования, реализованных в GROMACS.

Выбор силового поля для моделирования растворов литиевых солей на примере 1 М раствора LiClO₄ в сульфолане

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

Расчеты физико-химических свойств сульфолана, выполненные с применением различных силовых полей, показали, что ни одно из полей не позволяет одновременно рассчитать близкие к экспериментальным значения диэлектрической проницаемости, плотности и вязкости (табл. 1). Рассчитанное значение диэлектрической проницаемости сульфолана наиболее близко к экспериментально определенному при применении силовых полей Amber и Opslaa, плотности – Oplsaa, а вязкости – Gromos и Charmm.

В табл. 2 суммированы усредненные данные рассчитанных значений физико-химических свойств 1 М раствора перхлората лития в сульфолане c применением силовых полей, входящих в пакет GROMACS. Наиболее близки к измеренным рассчитанные величины плотности и вязкости 1 М раствора LiClO₄ в сульфолане при использовании силовых полей Gromos и Amber (табл. 2). При использовании силовых полей Opslaa и Charm рассчитанные величины плотности и вязкости отличаются не более, чем на 30% от измеренных. Также близки к измеренным (меньше в 2–3 раза) значения электропроводности и коэффициентов диффузии, рассчитанные при использовании силового поля Opslaa. Таким образом, наилучшие результаты моделирования физико-химических свойств 1 М раствора LiClO₄ в сульфолане методом МД достигаются при применении силового поля Opslaa.

Структура сольватных оболочек катиона лития

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

Выбор силового поля при моделировании 1 М раствора перхлората лития в сульфолане оказывает существенное влияние на форму ФРР (рис. 2а). Так, например, положение первого пика (и его высота) на кривой ФРР, описывающей распределение атомов кислорода молекул сульфолана вокруг катионов лития, изменяются следующим образом: Amber – 0.198 нм (45), Opslaa – 0.192 нм (42), Charmm – 0.146 нм (15) и Gromos – 0.414 нм (2.8).

Рис. 2.

Влияние силового поля на функции радиального распределения атомов кислорода, входящих в молекулы сульфолана (а) и аниона перхлората (б), вокруг катиона Li⁺ в 1M LiClO₄ в сульфолане; 1 – Amber, 2 – Opslaa, 3 – Charmm, 4 – Gromos.

Известно, что первая сольватная оболочка катиона лития обычно располагается на расстоянии 1.7–2.1 нм [2]. Расчетные величины расстояния между катионом лития и атомами кислорода сульфолана при использовании силовых полей Amber и Opslaa соответствуют этому диапазону. Расстояние между катионом лития и атомами кислорода молекул сульфолана, рассчитанное с применением силового поля Gromos, составляет 0.41 нм, что существенно больше расстояния между катионом лития и его первичной сольватной оболочкой. В случае применения силового поля Charmm рассчитанное расстояние между катионом Li⁺ и атомами кислорода сульфолана меньше расстояния до первой сольватной оболочки катиона лития и составляет 0.146 нм.

Выбор силового поля также оказывает влияние на форму кривых ФРР, описывающих распределение атомов кислорода, входящих в состав аниона ${\text{ClO}}_{4}^{ - }$, вокруг катиона лития Li⁺ (рис. 2б). Положение первого пика на кривой ФРР Li– ${\text{ClO}}_{4}^{ - }$, рассчитанной на основе результатов моделирования МД в силовых полях Gromos и Charmm, существенно отличается от положения аналогичных пиков, рассчитанных в силовых полях Opslaa и Amber, и выходит за диапазон размещения первой сольватной оболочки катиона лития. Положения первых пиков и их высоты кривых ФРР Li–О(${\text{ClO}}_{4}^{ - }$) составляют при использовании силового поля Amber – 0.202 нм (15), Opslaa – 0.238 нм (2.7), Charmm – 0.144 нм (7.7) и Gromos – 0.412 нм (7.5).

Существенное отличие форм и положений максимумов кривых функций радиального распределения Li–О(SL) и Li–О(${\text{ClO}}_{4}^{ - }$), рассчитанных на основе результатов моделирования МД в силовых полях Gromos и Сharm, указывает на некорректное описание взаимодействий катионов лития c анионами перхлората и молекулами сульфолана.

На основе анализа рассчитанных величин физико-химических свойств сульфолана и 1 М раствора LiClO₄ в сульфолане, а также функций радиального распределения вокруг катиона лития атомов кислорода сульфолана и аниона перхлората для дальнейшей работы нами было выбрано силовое поле Opslaa.

Влияние распределения зарядов на атомных центрах компонентов раствора на физико-химические свойства 1 М раствора LiClO₄ в сульфолане

Характер распределения зарядов на атомных центрах растворителей и анионов неизбежно будет влиять на рассчитанные значения физико-химических свойств электролитных систем. Поэтому нами было изучено влияние величин парциальных зарядов атомных центров молекул сульфолана и анионов перхлората на результаты расчетов физико-химических свойств 1M раствора LiClO₄ в сульфолане. Использованные в расчетах величины зарядов атомных центров перхлоратного аниона и сульфолана приведены в табл. 3 и 4.

Следует отметить, что свойства сульфолана, рассчитанные в работах [3, 29, 31] с применением зарядов атомных центров, представленных в табл. 4, очень хорошо совпадали с измеренными величинами. Однако наши исследования показали, что при моделировании свойств 1 М раствора LiClO₄ в сульфолане значения зарядов атомных центров сульфолана оказывают более сильное влияние на результаты расчетов (табл. 5). Независимо от используемых величин зарядов атомных центров расчетные величины плотности 1 М раствора LiClO₄ в сульфолане близки между собой (разброс значений не превышает 4%) и близки к экспериментально измеренным (табл. 5). Однако при расчете транспортных свойств (вязкости, электропроводности и коэффициентов самодиффузии) разброс рассчитанных значений достигает 350–400%.

Наиболее близки к измеренным величинам физико-химических свойств 1 М раствора LiClO₄ в сульфолане (табл. 5) рассчитанные с зарядами на атомных центрах перхлората (представленными в колонке 2 табл. 3) и сульфолана (колонка 4 табл. 4). Для дальнейших расчетов использовали эти заряды на атомных центрах сульфолана и аниона перхлората.

В неполяризуемых силовых полях переоцениваются межионные взаимодействия [38–40]. Это приводит к занижению расчетных величин электропроводности и коэффициентов самодиффузии. Для нивелирования данного эффекта заряды катионов и анионов растворенных солей уменьшают обычно на 10–30% [38–40]. С целью более точного описания свойств и строения 1 М раствора перхлората лития в сульфолане нами было выполнено моделирование растворов методом МД с применением указанного подхода – заряды катиона лития и аниона перхлората были уменьшены до 70, 80 и 90%.

Уменьшение зарядов катиона и аниона привело к увеличению расчетных величин электропроводности и коэффициентов самодиффузии (табл. 6). Рассчитанная электропроводность 1 М раствора LiClO₄ в сульфолане и коэффициент самодиффузии катиона лития стали сопоставимы с измеренными величинами. Однако рассчитанный коэффициент самодиффузии сульфолана стал в 4–5 раз выше измеренного.

Электропроводность растворов перхлората лития в сульфолане, рассчитанная по уравнению Эйнштейна–Гельфанда, ниже рассчитанной по уравнению Нернста–Эйнштейна. Это может быть объяснено тем, что уравнение Эйнштейна–Гельфанда учитывает все взаимодействия: анион–анион, анион–катион и катион–катион, в то время как уравнение Нернста–Эйнштейна учитывает только самодиффузию ионов [37].

Независимо от доли учета зарядов катиона лития и аниона перхлората коэффициенты самодиффузии компонентов 1 М раствора LiClO₄ в сульфолане уменьшаются в ряду D_SL > ${{D}_{{{\text{Cl}}{{{\text{O}}}_{{\text{4}}}}}}}$ > D_Li.

Согласно уравнению Стокса–Эйнштейна, коэффициент диффузии частицы обратно пропорционален ее радиусу и вязкости раствора. Диаметр сольватного катиона лития больше диаметра сульфолана и аниона перхлората. Таким образом, тот факт, что коэффициент самодиффузии катиона лития меньше коэффициента диффузии сульфолана, указывает на то, что катион лития в 1 М растворе перхлората лития в сульфолане переносится в сольватированном виде.

Строение 1 M раствора перхлората лития в сульфолане (структура сольватных оболочек катиона лития)

Для более детального исследования строения раствора перхлората лития в сульфолане были построены функции радиального распределения атомов кислорода и серы сульфолана и атомов кислорода и хлора перхлоратного аниона вокруг катиона лития.

Учет заряда катиона лития и аниона перхлората не оказывает влияния на форму функций радиального распределения атомов кислорода и серы сульфолана вокруг катиона лития (рис. 3), поэтому ФРР показаны только для 100%-ного учета заряда. На ФРР атомов кислорода сульфолана вокруг катиона лития наблюдаются два пика (рис. 3): первый узкий четкий при 0.19 нм и второй уширенный при 0.41 нм. Расстояние между положениями максимумов пиков (0.22 нм) сопоставимо с расстоянием между атомами кислорода в молекуле сульфолана (0.256 нм). Это указывает на то, что второй нечеткий уширенный пик на кривой ФРР Li–O(SL) обусловлен вторым атомом кислорода той же молекулы сульфолана. Поскольку молекула сульфолана может ориентироваться под разными углами по отношению к катиону лития, второй пик – уширенный с нечетким максимумом.

Рис. 3.

Функция радиального распределения атомов O и S сульфолана (а) и ее интеграл (б).

На кривой ФРР Li–S(SL) наблюдается только один узкий максимум при 0.33 нм, положение которого смещено относительно положения максимума ФРР Li–О на 0.14 нм, что сопоставимо с длиной связи S–O (0.148 нм). Это подтверждает то, что катион лития взаимодействует лишь с одним атомом кислорода сульфолана, а второй атом кислорода не участвует в сольватации катиона лития и обращен в противоположную от него сторону.

ФРР второго атома кислорода сульфолана вокруг катиона лития полностью совпадает с ФРР первого атома кислорода сульфолана (рис. 3). Это указывает на одинаковую вероятность сольватации катиона лития первым или вторым атомом кислорода сульфолана, что подтверждается равнозначностью рассчитанных координационных чисел катиона лития по атомам кислорода (табл. 7). Координационное число катиона лития по каждому атому кислорода составило 2.67 при 100%-ном учете зарядов ионов (табл. 7). Общее координационное число катиона лития по обоим атомам кислорода составляет 5.34, а координационное число Li⁺ по атомам серы – 5.35. Таким образом, наиболее вероятное расстояние Li–О в районе 0.2 нм совместно с координационным числом, равным 5.3, свидетельствует об образовании достаточно структурированной первой сольватной оболочки катиона Li⁺ в растворе перхлората в сульфолане. Рассчитанные величины координационного числа лития согласуются с литературными данными. Обычно координационное число катиона лития в электролите изменяется в диапазоне 4–5 [41]. Уменьшение доли заряда ионов привело к снижению рассчитанных координационных чисел катиона лития по атомам кислорода и серы сульфолана (табл. 7).

Формы ФРР различных атомов кислорода перхлоратного аниона различны (рис. 4): отличаются высоты пиков, максимумы которых расположены на расстояниях 0.21 и 0.46 нм от Li⁺, а высота третьего пика ФРР (при 0.76 нм) примерно одинакова для всех атомов кислорода перхлоратного аниона. На кривой ФРР атома хлора вокруг катиона лития находятся два пика, максимумы которых расположены при 0.36 и 0.66 нм. Это может указывать на то, что перхлоратный анион может входить как в первую, так и во вторую сольватные оболочки катиона лития.

Рис. 4.

Радиальные функции распределения атомов кислорода и хлора перхлоратного аниона вокруг катиона лития в 1 М растворе LiClO₄ в сульфолане (а, в) и их интегралы (б, г).

Учет зарядов ионов оказывает существенное влияние на формы ФРР атомов кислорода и хлора аниона перхлората вокруг катиона лития (рис. 4). Соотношение высот первого пика с высотами второго и третьего на кривых ФРР O(${\text{ClO}}_{4}^{ - }$)–Li уменьшается с уменьшением доли заряда ионов в растворе. В случае 70 и 80% от зарядов ионов на кривой функции радиального распределения первый и второй пики практически исчезают. Аналогичным образом исчезает первый пик на кривой ФРР Cl(${\text{ClO}}_{4}^{ - }$)–Li. С уменьшением доли заряда ионов взаимодействие катиона лития и аниона перхлората лития ослабевает.

На кривых интеграла функций радиального распределения O(${\text{ClO}}_{4}^{ - }$)–Li отсутствуют плато (рис. 4). Это указывает на то, что расположение перхлората в первой сольватной оболочке катиона лития не структурировано в отличие от молекул сульфолана.

Уменьшение доли зарядов ионов привело к снижению рассчитанных координационных чисел катиона лития по атомам кислорода и хлора аниона перхлората (табл. 8).