Open: Мастер-класс: Моделирование реакции

Введение

В данном мастер-классе/мануале я покажу вам как смоделировать реакцию от экспериментальных данных до нахождения механизма реакции.

Что будет в мастер-классе:

• Составление ментальной модели реакции
• Предложение механизма реакции
• Выбор приближений
• Построение структур в Chemcraft и составление инпут-файлов для ORCA
• Запуск расчетов и анализ получившихся файлов
• (Возможно) Нестандартные типы расчетов (IRC, NEB)

Чего не будет

• Конформационного поиска (по нему есть свой мануал, и не один)
• Перевода с теоретического на химический, анализ взаимодействий и орбиталей (тоже есть свой мануал)

Реакция взята отсюда [pdf], все данные в статье есть, поэтому если захотите найти подсказки, всегда сможете это сделать. В самом мастер-классе я постараюсь свести подсказки к минимуму и даже в конце не приведу правильный механизм, посмотрите его в статье

Даже если вы не будете читать дальше, скажу самое главное в моделирование реакций:

• Внимательно изучите экспериментальные данные
• Читайте мануал используемой программы
• Иногда экспериментальные данные могут ошибочны, а экспериментаторы могут неправильно их интерпретировать
• Иногда расчеты могут показывать неверные данные, а расчетчики могут неправильно их интерпретировать

Что понадобится

• Доступ к боту/серверву
• Chemcraft
• Notepad++
• Chemdraw
• (Опционально) Python/Excel/Photoshop для построения энергетического профиля реакции

Основы моделирования

По нашему мнению любое моделирование состоит из следующих этапов (это все изложено в лекции, тут)

1. Построение ментальной модели процесса:
   • Вы должны добыть и осознать всю доступную информацию про реакцию. Вообще все что могут вам сказать экспериментаторы. перечислю самые основные данные
     • Условия эксперимента (кто к чему доливали, в каком порядке, какие условия, ...), а так же всю информацию про оптимизацию условий (перебор р-лей, температур, ...)
     • Наблюдения в процессе реакции (изменения цвета, выделение газа, нагрев смеси, ...) и после нее (выход, конверсия, чистота, ...)
     • Данные ф/х методов анализа (ЯМР, ИК, рентген) и в какой момент они сняты (в процессе реакции, через час/два/день)
     • Известные химикам релевантные статьи по теме
     • Мнение химиков обо всем этом. При этом помните, что ваши коллеги это тоже люди и ученые, которые могут ошибаться. Относитесь критично к их анализу
   • Теперь начинается ваша работа. И начинается она с изучения литературы. Поищите сами статьи по теме проекта (и экспериментальные и расчётные)
   • После этого можно наконец проанализировать все данные и
     • Понять на какие вопросы будет отвечать ваш фрагмент статьи (Например, вы объясняете ee или необычный продукт)
     • Предложить механизм процесса и какие эффекты его контролируют. Обязательно нарисовать механизмы в chemdraw. Да именно множественное число, сразу разрисуйте все ваши теории про протекание реакции.
2. Перенос ментальной модели в расчет. Поздравляю вы предположили что происходит в колбе. Теперь надо понять как это посчитать. Напомню, что у нас есть 3 основных класса приближений. Будет кратко, подробнее есть в лекции. Еще раз скажу, что на многие вопросы можно найти ответы в статьях. С большой вероятностью подобные реакции кто-то уже моделировал. Но опять-таки, относитесь критично к написанному: если в статье приведено сравнение 3 методов для задачи и обосновано выбран один этому в целом можно верить, а если люди просто взяли B3LYP/6-31G, то подумайте дважды.
   • Химическая модель (какие атомы во вселенной учитываем и как).
   • Ансамбль (конф. поиск, молекулярная динамика, ...)
   • Уровень теории (как решаем уравнение Шредингера)
3. Рисование структур в Chemcraft (или не дай бог где-то еще) и постановка расчетов
4. Анализ результатов и либо написание статьи, либо возвращение к пункту 1

Теперь давайте поговорим подробнее про пункты 3 и 4

Chemcraft

Большая часть популярных программ для квантовой химии использует формат координат xyz. То есть трехмерные координаты в виде:

Символ элемента (или порядковый номер) координата по x координата по y координата по z

C       -2.727421687      0.000000000     -5.280900664
H       -2.196987741     -0.951333822     -5.204710633
H       -2.010874002      0.821859802     -5.223394016
H       -3.443969372      0.083844265     -4.461308883
H       -3.257855632      0.045629756     -6.234189124

Сложность заключается в том, что есть по факту два формата xyz: обрезанный (в примере выше) и полный (ниже), в котором первая строка это число атомов, вторая комментарий (любой или пустая строка) и только затем координаты. Некоторые программы (orca, gaussian) требуют обрезанный формат, а другие полный (xtb, crest). Следите за этим

5
This is CH4 in full xyz
C       -2.727421687      0.000000000     -5.280900664
H       -2.196987741     -0.951333822     -5.204710633
H       -2.010874002      0.821859802     -5.223394016
H       -3.443969372      0.083844265     -4.461308883
H       -3.257855632      0.045629756     -6.234189124

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

Номера атомов в Chemcraft начинаются с "1", как в большинстве программ, но в ORCA нумерация идет с "0", поэтому, например, при постановке констрейнов (дальше рассказано про это) в инпут файле для ORCA, вычитайте единицу из номера атома в Chemcraft

Учить пользоваться Chemcraft через текст задача неблагодарная. Поэтому я ее выполнять не буду, но дам пару советов.

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

В структуре должны быть правильные длины связей и углы между ними, большая часть связей должна находиться в заторможенной конформации (если это возможно). Следите за расстояниями между несвязанными атомами. Например, может "случайно" образоваться молекула H₂ из близко расположенных атомов метильных групп или другие подобные недоразумения. Важно запомнить еще три вещи про рисование структур.

• Химическая связь это модель, и в реальности она не существует, большая часть квантовохимических программ оперирует только длинами связей, если вы не указываете обратного и не задаёте порядок связи
• Программа не способна изменять конфигурацию E/Z изомеров и находить оптимальные торсионные углы. Если в вы не можете опрделить их вручную, используйте конформацонный поиск (CREST, Maestro)
• При открытии результата расчета в Chemcraft могут не отображаться некоторые связи из-за их большой длины (зачастую связи с металлом, и еще с эпоксидами у него проблемы), в таком случае посмотрите на длины связей и оцените их адекватность
  • При этом Chemcraft все еще очень умная программа и 95% связей он рисует правильно, этим можно воспользоваться: нарисуйте вашу структуры как считаете правильным, сохраните файл, закройте Chemcraft, а затем откройте структуру заново. Chemcraft самостоятельно проставит везде связи и если вы увидите, что он поставил связь там, где вы их не ждете (пример про H₂ выше) или не нарисует связь там, где ждете, то с очень большой вероятностью вы нарисовали структуру криво, подумайте и перерисуйте.

Чтобы нарисовать наиболее правдоподобную геометрию можно изучить структуры из РСА похожих соединений и статьи про них (просите у химиков либо в CCDC). Помните однако, что структура в кристалле не всегда соответствует структуре в растворе.

Для начала я всегда советую вам сохранять нарисованную структуры в формате xyz, пускай лежат.

Задача

Изучаемая реакция

Есть два соединения, содержащих хиназолинононовый и цимантренильный фрагменты (4 и 3). В ацетонитриле происходит отщепление одного CO-лиганда и реакция получающегося шестнадцатиэлектронного комплекса с молекулой растворителя. В бензоле соединение 4 реагируют в соответствии с нашими ожиданиями, в результате чего образуется хелат 6, со связью Mn-O. В то же время реакция с соединением 3 не приводит к продукту 5, а вместо этого получается смесь 4 и 6. Это подтверждается спектрами ИК, ЯМР и УФ. Насколько нам известно, этот фотохимический перенос цимантренильной группы от кислорода к азоту является первым примером, найденным в литературе, где цимантрен сохраняет все три карбонила при фотовозбуждении. Надо понять, почему это так происходит, то есть предложить механизм реакции из 3 в 4. Предложить механизм, значит локализовать все минимумы (в расчете не должно быть мнимых частот) и переходные состояния (ПС) (1 мнимая частота, соответствующая разрыву и образования связей в ходе реакции), а также найти все энергии активации (показать, что константы скорости соответствуют эксперименту) и свободные энергии Гиббса (показать, что процесс разрешен).

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

1. Процесс термической природы с миграцией CH₂ от азота к кислороду
2. Реакция с сохранением всех трех CO-лигандов и участием марганца в реакции
3. Реакция, начинающаяся с удаления одного CO-лиганда с перегруппировкой и возвращением CO

Порядок действий

Вам нужно нарисовать в Chemdraw предполагаемый механизм, по одной (или нескольким) моделям. Он должен начинаться с 3 и заканчиваться 4. Постарайтесь применить свои знания химии и посмотреть в статьях вероятные структуры интермедиатов. Я прикрепляю файл со своими механизмами.

Далее рисуем структуры в Chemcraft

Теперь надо выбрать приближения

• Метод решения уравнения Шредингера (не забыть про полуэмпирическую поправку)
• Базис
• Метод учета сольватации
• Метод расчета энергии (свободная/электронная)

Если это сделано, можно создать инпут-файл. Прикрепляю пример, который получился у меня, с ним можно сверить ваш. Оптимизация в минимум, Поиск переходного состояния, Оптимизация с констрейном (не забывайте делать ее перед поиском ПС, а также помните, что в ORCA нумерация атомов начинается с нуля)

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

!! Далее я буду приводить советы по расчетам, но не давать подробное описание упоминаемых keywords, читайте мануал или спрашивайте !!

1. Расчет завершился успешно, геометрия выглядит, как мы ожидаем, в минимуме нет мнимых частот, в переходном ― одна мнимая частота, отвечающая реакции. В таком случае можно смотреть на энергии (Final Gibbs free energy), найти энергии активации и оценить константы скорости реакции, чтобы сделать вывод о возможности протекания реакции через исследуемый механизм. Напомню, что при комнатной температуре энергия активации не должна превышать 25 ккал/моль, чтобы реакция шла. Есть четыре возможных варианта, при которых энергия активации на одной из стадий получается больше 25 ккал/моль.
   • Слишком высокая энергия максимума (ПС). Основная причина - не минимальная по энергии конформация
   • Слишком низкая энергия минимума. Вы спросите, может ли такое быть? Да, может. Структура, которую вы нашли как минимальную может не существовать в реальности. Например, вы координировали лиганд на металл бидентатно, хотя в реальности одну из вакансий занимает растворитель, что вы можете проверить рассчитав комплекс с растворителем.
   • Неправильно рассчитаны энергии. То есть подвело одно из приближений. Может базис слишком маленький или нужен явный учет растворителя.
   • Вы посчитали все правильно, но модель реакции нерелевантная и не может описать поведение системы в этих условиях.
2. Расчет упал и выдал техническую ошибку. Ошибок в расчете может быть великое множество и обсуждать их все я не хочу. Если вы не понимаете смысл ошибки, можете мне писать, да и в целом, если хотите моделировать реакции присоединяйтесь к друзьями нашей группы, я устраиваю созвоны раз в неделю, отвечая на вопросы. Частыми причинами ошибок могут стать: плохая начальная геометрия, неправильный заряд и мультиплетность, ошибка в инпуте, недостаток памяти (можно увеличить параметром %maxcore 3000, то есть выделить на каждое ядро по 3ГБ оперативной памяти). Исправив ошибку переставьте расчет
3. Расчет свелся куда-то не туда или вообще не свелся. Тут опять-таки может быть множество проблем, поэтому дам лишь пару советов
   • Плохая начальная геометрия. Вы можете сделать следующее
     • Взять последнюю геометрию из расчета, если она похожа на ожидаемую
     • Сделать предоптимизацию
     • Перестроить исходя из советов выше
   • Вы уже находитесь очень близко к минимуму или максимуму, и/или критическая точка не выражена. То есть расчет постоянно "перепрыгивает" правильную геометрию. В таком случае имеет смысл поварьировать параметр %maxstep, изменяя шаг расчета.
   • Часто при поиске ПС помогает пересчет матрицы Гесса через каждые несколько шагов, эту функцию можно включить параметром Recalc_Hess
   • Того, что вы ищите, не существует. Программа может не найти минимум или ПС, потому что его не существует даже на бумаге. Как отличить этот исход от предыдущих? Короткий ответ - простого метода не существует. Для этого нужен опыт в теоретической химии и знания в синтетической. Кроме того обычно приходится перепробовать несколько из вышеописанных (и неописанных) методов сведения расчета, чтобы более-менее убедить себя, что искомой структуры нет. Самый "простой" способ выйти из этого положения - найти правильную структуру, согласующуюся со всеми данными или полностью опровергнуть модель, на которой построен поиск этого несуществующего соединения
   • Попробуйте воспользоваться другим методом поиска ПС или минимума. Во-первых, примените оптимизацию с констрейном, если не сделали этого. Если ищите ПС и у вас есть реагент и продукт можно сделать NEB (Nudged Elastic Band). Этот метод "соединяет" два минимума, пытаясь найти максимум на потенциальной кривой, тем самым получая ПС. Если у вас получилось какое-то ПС, но мнимая частота не совсем "понятная", то есть не очевидно, какие два минимума соединяет это ПС можно сделать IRC (Intrinsic Reaction Coordinate), которое как раз находит реагент и продукт, соединённые вашим ПС

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

Профиль реакции

Еще немного советов:

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

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

Спасибо мне (Чалому Василию) за то, что я это делаю.

Александре Дворжанской за редакторскую правку этого мастер-класса. Благодаря ей этот текст можно прочитать и понять.