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

Введение

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

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

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

Чего не будет

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

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

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

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

Задача

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

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

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

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

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

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

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

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

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

• Химическая связь это модель, и в реальности она не существует, большая часть квантовохимических программ оперирует только длинами связей, если вы не указываете обратного и не задаёте порядок связи
• Программа не способна изменять конфигурацию E/Z изомеров и находить оптимальные торсионные углы. Если в вы не можете опрделить их вручную, используйте конформацонный поиск (CREST, Maestro)
• При открытии результата расчета в 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), которое как раз находит реагент и продукт, соединённые вашим ПС

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

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

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

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

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

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