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

ЭКСПЕРИМЕНТАЛЬНАЯ ЧАСТЬ

Методы

В работе рассматриваются таутомерные формы нейтральной молекулы 4,5-диметил-2-(2'-гидроксифенил)имидазола в водном растворе. Согласно литературным данным, в воде при pH 8.4 присутствуют три таутомерных формы: цис-енольная (E_c), транс-енольная (E_t) и кето- (K) формы (рис. 1) [13]. При построении модельной системы первоначально проводилась оптимизация геометрических параметров ДМГИ методом функционала электронной плотности в варианте PBE0/6-31G**. Этот протокол хорошо описывает равновесные геометрические параметры, а также относительные энергии конформеров молекул [6], в том числе органических хромофоров, стабилизированных водородными связями [14]. Расчеты всех систем проводились в основном синглетном состоянии, S₀.

Рис. 1.

Таутомерные формы нейтральной молекулы 4,5-диметил-2-(2'-гидроксифенил)имидазола: цис-енольная (E_c), транс-енольная (E_t) и кето- (K).

Далее молекулу ДМГИ помещали в прямоугольный параллелепипед из молекул воды так, чтобы расстояние от органической молекулы до границ ячейки составляло не менее 10 Å. Размер модельной системы составил 1394 атома в ячейке размером 30 × 30 × 30 ų (рис. 2).

Рис. 2.

Модельная система для расчетов методом КМ/ММ. Здесь и далее шаростержневое представление относится к ДМГИ, стержнями показаны молекулы воды в КМ-подсистеме, линиями – молекулы воды в ММ-подсистеме. Атомы азота и кислорода, участвующие в таутомеризации маркированы O1, N1 и N2. Пунктирной окружностью выделен атом, от которого рассчитывается расстояние до молекул воды при включении их в КМ- или ММ-подсистему.

Для первичного уравновешивания модельной системы проводились расчеты классической молекулярно-динамической траектории длиной 50 пс. Классические и КМ/ММ МД-расчеты проводились в NPT-ансамбле, при температуре 300 К и давлении 1 атм с использованием термостата Ланжевена с параметром связывания 50 фс^–1 и баростата Нозе–Гувера [15], шаг интегрирования составлял 0.5 фс. Взаимодействия рассчитывались для атомов, находящихся на расстоянии менее чем 12 Å в классическом расчете и для ММ-подсистемы в КМ/ММ расчете. В ходе классической молекулярной динамики ДМГИ описывался набором параметров силового поля CGenFF [16], а молекулы воды – параметрами TIP3P [17]. Эти расчеты проводились в программном пакете NAMD 2.12 [12]. Классические силовые поля не позволяют описывать таутомерные переходы, то есть акты образования и разрыва химических связей, поэтому все последующие расчеты проводились методом молекулярной динамики с потенциалами КМ/ММ. Здесь стоит отметить особенность рассматриваемых систем: молекулы воды очень лабильные и в ходе расчета траектории могут входить в ближайшее окружение ДМГИ и выходить из него. Иными словами, если в начале моделирования жестко фиксировать атомы, которые будут входить в КМ-подсистему, то в ходе молекулярно-динамического расчета КМ- и ММ-молекулы воды перемешаются, что повлияет на описание взаимодействий между растворенной молекулой и растворителем. Поэтому в данной работе за основу был взят алгоритм, разработанный авторами пакета NAMD для проведения КМ/ММ МД-расчетов, позволяющий переопределять КМ-подсистему в ходе расчетов, руководствуясь определенными правилами [18]. Интерфейс, предложенный разработчиками, предполагал взаимодействие квантово-механического пакета ORCA [19] c МД-пакетом NAMD [12]. В нашей работе мы создали новый интерфейс, позволяющий проводить расчеты методом КМ/ММ МД, используя для расчета энергий и сил в КМ-подсистеме программный пакет Firefly [10]. Это открывает ряд новых возможностей, в частности, эффективные высокопараллельные расчеты на суперкомпьютерах, а также применение всего арсенала многоконфигурационных методов и наиболее точных на сегодняшний день вариантов теории возмущения, реализованных в Firefly [20].

При проведении КМ/ММ МД-расчетов КМ-подсистема включала в себя молекулу ДМГИ и молекулы воды, у которых хотя бы один из атомов оказывался на расстоянии 7 Å от атома углерода, находящегося между атомами азота в имидазольном кольце (рис. 2). Такое условие позволило включать только ближайшие к ДМГИ молекулы воды, которые образуют водородные связи с атомами азота и кислорода, участвующими в процессе таутомеризации, а также формировать цепочку водородных связей для переноса протона между дистанцированными атомами. Комбинация NAMD/Firefly позволяет на каждом шаге МД-траектории обновлять состав молекул растворителя в КМ-подсистеме, поэтому состав растворителя менялся при переходе от кадра к кадру в зависимости от расстояния до растворенной молекулы.

Молекулярно-динамические расчеты начинались с трех модельных систем, находящихся в различных таутомерных формах, длина каждой из траекторий составила 4 пс. Расчеты проводились методом PBE0/6-31G**/TIP3P. Дополнительно получена молекулярно-динамическая траектория большей длины в минимальном базисе 3‑21G. Это было необходимо для наблюдения заформированием сетки водородных связей между атомами N1 и O1, расположенными с противоположных сторон хромофора. Выбранный вариант КМ/ММ позволяет учитывать электростатические взаимодействия КМ- и ММ-подсистем. Для КМ-подсистемы вклады от ММ-атомов учитываются за счет вкладов от точечных зарядов, предписанных силовым полем, в одноэлекторнную часть гамильтониана квантовой подсистемы. Нековалентные взаимодействия на ММ-уровне учитываются за счет электростатического члена, рассчитывающегося из зарядов ММ-атомов и зарядов Малликена атомов КМ-подсистемы, а также ван-дер-ваальсовых взаимодействий, параметры для которых для всех атомов системы берутся из силового поля [16].

Для описания спектров поглощения всех таутомерных форм использовались молекулярные кластеры, полученные путем выделения из различных кадров молекулярно-динамической траектории молекулы ДМГИ вместе с первой сольватной оболочкой. Энергии вертикальных переходов в первое возбужденное синглетное состояние рассчитывались нестационарным методом функционала электронной плотности (TDDFT) c функционалом ωB97x-D3 в базисе def2-SVP в программном пакете ORCA [19]. Такой подход удовлетворительно показал себя при описании S₀–S₁-переходов для органических хромофоров [21–23]. Известно, что метод TDDFT дает завышенные энергии электронных переходов, поэтому обычной практикой является добавление аддитивного члена для лучшего согласия расчетов и эксперимента [24]. В нашем случае мы вычитали 1000 см^–1 из всех получаемых в расчете значений. Для каждой из таутомерных форм расчеты энергий вертикальных переходов проводились в 50 точках траектории.

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

Расчеты методом молекулярной динамики с КМ/ММ потенциалами позволяют наблюдать акты образования и разрыва химических связей, в нашем случае таутомерные превращения. На рис. 3 показан пример переноса протона от атома N2 на атом O1 в ходе расчета (здесь и далее названия атомов соответствуют рис. 3). В качестве исходной структуры использовалась кето-форма, K, протонированная по двум атомам азота. В ходе молекулярно-динамического расчета протон переходит на атом кислорода O1 с образованием цис-енольной формы, E_c.

Рис. 3.

Распределение расстояний O1···H (серым) и N2···Н (черным) вдоль молекулярно-динамической траектории, демонстрирующее таутомерный переход между цис-енольной и кето-формами в ходе молекулярно-динамического моделирования.

Важно отметить значимость учета водного окружения при моделировании системы. Согласно экспериментальным данным по спектрам поглощения [13], которые будут подтверждены расчетами энергий вертикальных электронных переходов далее в тексте, кето-форма преобладает в равновесии над нейтральными формами ДМГИ. Тестовые расчеты ДМГИ в газовой фазе показали невозможность локализовать минимум на поверхности потенциальной энергии, соответствующий кето-форме K. Добавление двух молекул воды, образующих водородные связи с атомами O1 и N1 органической молекулы, также не позволили локализовать минимум, соответствующий форме K.

Другое обстоятельство связано с корректной оценкой относительных энергий конформеров. Согласно нашим расчетам в газовой фазе, форма E_t на 7 ккал/моль выше по энергии, чем E_c, что соответствует константе равновесия порядка 10^–5 при комнатной температуре и 10^–16 при 100 K. Добавление только двух молекул воды снижает разность энергий до 3 ккал/моль, что уже соответствует более реалистичным значениям. Дальнейшее расширение модельной системы приводит к появлению множества конформаций молекул растворителя, и становится невозможным сравнивать относительные энергии таутомеров. В таких случаях принято использовать метод молекулярной динамики с ограничивающим потенциалом, позволяющий проводить серию расчетов, в которых сканируются области конформационного пространства вблизи заданного значения координаты реакции [25]. Эти расчеты являются крайне ресурсоемкими, поэтому в данной работе мы изучали принципиальную возможность образования сетки водородных связей для переноса протона между атомами N1 и O1. Для этого мы проводили более продолжительный расчет, длина траектории составила 40 пс, и использовали менее затратный метод описания КМ-подсистемы – PBE0/3-21G. На рис. 4 показана цепочка из молекул воды, позволяющая осуществить перенос протона с атома N1 на атом O1. Она образуется через 15 пс после начала расчета и далее остается в основном (около 80% времени) стабильной в ходе траектории. Эти результаты могут быть использованы в качестве стартового приближения для поиска координаты реакции переноса протона.

Рис. 4.

Сетка водородных связей, позволяющая осуществить перенос протона между атомами N1 и O1.

В работе [13] представлены данные электронных спектров для молекулы ДМГИ в водном растворе при pH 8.4, однако соотнесение экспериментальных полос поглощения с таутомерными формами без проведения расчетов является затруднительным. Мы использовали три молекулярно-динамические траектории с различными таутомерами ДМГИ и случайным образом выбирали по 50 кадров. Системы редуцировались так, чтобы вокруг растворенной молекулы была одна сольватная оболочка, и проводили расчет энергии вертикального S₀–S₁-перехода. Далее проводилась статистическая обработка результатов. Полученные значения кластеризовались с шагом 500 см^–1, после чего строились распределения по волновым числам. После чего распределения смещались на 1000 см^–1, чтобы добиться лучшего согласия с экспериментом. На рис. 5 представлено наложение экспериментального спектра, полученного оцифровкой спектра поглощения в водном растворе при pH 8.4 на рис. 1б в работе [13] и рассчитанных спектров для каждой из таутомерных форм. Обе енольные формы характеризуются полосами поглощения в области 33 000–39 000 см^–1. Для кето-формы характерны меньшие значения волновых чисел в области 27 000–32 000 см^–1. Наложение рассчитанных и экспериментальных спектров позволяет провести отнесение полос соответствующим таутомерным формам. Максимум полосы поглощения около 35 000 см^–1 относится к цис-енольной форме, а полоса с максимумом в районе 36 500 см^–1 – к транс-енольной. Таким образом, два максимума широкой полосы в этой области соответствуют двум таутомерам. Согласно данным расчетов, широкую полосу с двумя максимумами в области меньших волновых чисел стоит отнести к поглощению кето-формы, а наличие двух максимумов можно объяснить двумя популяциями несколько отличающихся друг от друга конформаций.

Рис. 5.

Экспериментальный спектр поглощения ДМГИ в водном растворе при pH 8.4 (серая линия, 1) [13] и рассчитанные полосы поглощения для кето- (черная сплошная линия, 2), цис-енольной (черная штрихпунктирная линия, 3) и транс-енольной (черная пунктирная линия, 4) форм.

Чтобы объяснить относительные положения максимумов полос поглощения цис- и транс-енольных форм, мы провели дополнительные расчеты в газовой фазе. В рамках тех же вычислительных протоколов рассчитаны энергии вертикальных электронных переходов и разностные электронные плотности между основным и первым возбужденным электронным состояниями. По результатам расчетов энергия вертикального S_0,min–S₁-перехода составила 37 743 см^–1 для E_с-формы и 39 398 см^–1 для E_t-формы. Относительное расположение энергий вертикальных электронных переходов в газовой фазе воспроизводит таковое для рассчитанных полос поглощения в растворе (рис. 5). Проанализируем карты разностной электронной плотности (рис. 6). Для обеих систем характерно уменьшение электронной плотности на атоме кислорода O1 при переходе в S₁-состояние. В случае транс-конформера присутствует сильная внутримолекулярная водородная связь N2–H···O1, которая стабилизирует основное состояние и дестабилизирует возбужденное; в цис-конформере такой связи нет. Это и приводит к синему сдвигу E_t – относительно E_c-формы.

Рис. 6.

Разностная электронная плотность (S₁–S₀) для E_c- и E_t-форм. Плюсами показаны области увеличения электронной плотности при переходе из S₀- в S₁-состояние, минусами – области уменьшения.

Таким образом, применение метода молекулярной динамики с КМ/ММ потенциалами позволило описать переходы между таутомерными формами нейтральной молекулы ДМГИ в водном растворе. На основании рассчитанных спектров поглощения проведена интерпретация экспериментального спектра поглощения и проведено отнесение полос. Анализ разностных электронных плотностей при переходе из основного в первое возбужденное электронное состояние позволило объяснить синий сдвиг транс-енольной формы относительно цис-енольной формы ДМГИ.

Исследование выполнено при финансовой поддержке Российского фонда фундаментальных исследований в рамках научного проекта № 18-03-00605. Работа выполнена с использованием оборудования Центра коллективного пользования сверхвысокопроизводительными вычислительными ресурсами МГУ имени М.В. Ломоносова.