Методы компьютерного моделирования в химии

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

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

В молекулярной химии широко распространены методы расчётов и оптимизации геометрии молекул, расчетов их свойств, возможных колебаний и термодинамических параметров и др.

Оптимизация геометрии сложных молекул – один из самых широко используемых методов моделирования, который необходим для того, чтобы находить наиболее вероятные существующие конформации молекул, особенно в сложных белках и органической химии. Основан метод на всем известном со школы факте, что «всё стремится к минимуму свободной энергии», и сводится метод к тому, чтобы найти такие координаты атомов, при которых полная энергия молекулы минимальна. Вариаций и программных реализаций оптимизации геометрии огромное количество, а самыми популярными являются метод BFGS (алгоритм Бройдена-Флетчера-Голдфарба-Шанно) и метод сопряженных градиентов. Они основаны на градиентных вычислениях зависимостей полной энергии от изменения, так называемых, избыточных внутренних координатах атомов, т.е. межатомных расстояний, валентных и диэдральных углах. На входе алгоритм принимает начальную геометрию молекулы, которая представляет собой большой массив данных о расположении химических элементов, типах связи, их силе и т.п., а результат моделирования представляет собой новый массив, содержащий оптимизированную геометрию молекулы, полную энергию, производные энергии (градиенты), которые могут пригодится для дальнейших расчётов частот колебаний или энергии химической реакции. Для реализации на практике данных методов используется специализированное программное обеспечение, например, GUI Avogadro.

После оптимизации геометрии молекулы при наличии экспериментальных данных, можно произвести моделирование молекулярных свойств. Такими свойствами являются, например: орбитальные энергии, потенциалы ионизации, атомные заряды, дипольные и высшие мультипольные моменты молекулы; поверхности электростатического потенциала молекулы, электронной плотности молекулярных орбиталей, энергии электронно-возбужденных состояний, поправки к энергии молекулы, обусловленные влиянием среды и многое другое. Например, расчёт энергии электронных возбуждений производится с помощью теории функционала плотности (DTF), зависящего от времени. Этот метод из квантовой механики, который используется для численного решения большого количества приближённых уравнений Кона-Шэма. После его применения с помощью программ визуализации можно увидеть график спектра поглощения, т.е. зависимость поглощения фотонов от длинны волны.

Также одним из часто необходимых методов моделирования является нахождение колебательных частот молекулы. Они, в частности, определяют устойчивость молекулы, вид её поверхности потенциальной энергии (ППЭ), определяет точки переходных состояний в реакциях. Находят колебательные частоты с помощью составления матрицы Гесса, каждый элемент которой представляет собой вторую производную энергии по координатам атомов. После вся матрица диагонализируется, а собственные её значения делятся на массы атомов, что позволяет получить величины частот в квадрате. Стоит отметить, что этот метод моделирования обычно проводят после оптимизации геометрии во избежание ошибок расчёта. Также теория и опыт показывают, что реальные значения частот отличаются на 5-10% от рассчитанных (они часто завышены) за счёт того, что в реальности потенциал между атомами не соответствует модели гармонического, т.е. частоты являются ангармоническими. Похожим образом, решением секулярного уравнения, находятся нормальные моды колебаний, что способствует изучению колебательных спектров молекулы, т.е. зависимости интенсивности поглощения от частоты в разных диапазонах. Для практической реализации данных методов создано много программ и средств визуализации, например, Moltran и Gaussian.

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

Обычно при исследовании химических реакций исследователи используют уже имеющиеся данные о всех физических свойствах взаимодействующих веществ, либо моделируют самостоятельно, используя вышеперечисленные методы. На этапе моделирования химических реакций необходимо найти термодинамические функции и параметры, которые позволят узнать, например, константу равновесия реакции или её энергию, что очень важно практически во всех исследованиях. Термодинамические функции энтальпии H, энтропии S, внутренней энергии U, свободной энергии Гиббса G, свободной энергии Гельмгольца F и некоторые более специфичные функции могут вычисляются на основе квантово-химических данных с использованием статистической термодинамики, если до этого геометрия молекул была оптимизирована, были найдены энергетические спектры и мультиплетность молекул, а также были найдены колебательные частоты. Конкретные формулы могут отличаться от метода к методу, но в основе используются:

методы DFT и TD-DFT, наличие разнообразия функционалов которого позволяет использовать его для молекул практически любого вида;
методы связанных кластеров CCSD(T) – особенно часто применяется для ван-дер-ваальсовых молекул, хорошо вычисляют энтальпии реакций;
методы теории возмущений Мёллера-Плессета второго и четвёртого порядков широко применяются для малых молекул;
метод бракнеровских пар, который даёт хорошие оценки для энтальпий небольших молекул;
методы конфигурационного взаимодействия CIS и CISD, MCSCF и CAS SCF, CAS PT2, CAS MP2, CAS CCSD требуют очень больших вычислительных мощностей, но обеспечивают рекордную точность;
полуэмпирические методы PM3 и его разновидности PM6, PM7, AM1, ZINDO и ZINDO/S, которые могут иметь и малую погрешности, и малые вычислительные затраты по сравнению с другими методами, но хорошо моделируют свойства органических и неорганических соединений.

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

сначала подобрать вид функционала и базиса;
произвести расчёт геометрии и частот всех реагентов;
вычислить сумму полных энергий молекул и термические поправки;
расчет величины термодинамических функций.

Для более их точного вычисления можно использовать специальные методы, которые обычно называют высокоточными экстраполяционными методами определения энтальпий, например, G3, G4, CBS, W1BD. Но точность рассчитанной энтропии таким образом повысить нельзя, но к счастью, точность моделируемых величин статистически близка к измеренным экспериментально.

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

Моделирование переходных состояний и скоростей реакций

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

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

Если количество точек переходного состоя одна или две, то производится оптимизация геометрии, а в ином случае необходимо заново вручную изменить геометрию переходного комплекса. Эта оптимизация проходит с помощью компьютерной реализации метода следования по собственному вектору, основанному на векторной алгебре. В найденной точке производится расчёт частот, а после найденные частоты анализируются. Если мнимая частота единственна, то это ещё не доказывает, что найдено нужное переходное состояние, ведь, например, реакция может быть многостадийной. Для решения этой проблемы необходимо установить путь реакции через данное состояние, что делается с помощью процедуры построения внутренней координаты реакции (IRC, intrinsic reaction coordinate). Если все расчёты закончились успешно и путь реакции был построен, то можно приступить к определению скоростей стадий реакции с помощью стандартной теории переходного состояния (CTST). Для этого необходимо знать термодинамические функции реагентов и переходного состояния (для этого нужно будет промоделировать и его). Конкретные формулы для вероятностей реакции и скоростям зависят от количества реагентов, их концентрации и фаз реагентов. Их можно найти, как и другие формулы, в учебниках по квантовой химии и термодинамике.

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

метод самосогласованного реакционного поля (например, модель Онзагера) – он описывает перераспределение атомных зарядов в молекуле, т.е. изменение поляризации, за счёт дополнительного потенциала на границе с конденсированной фазой;
модель поляризуемого континуума:
- PCM использует полости размеров ван-дер-ваальсовых радиусов, которая лучше описывает контакты молекулы и окружения;
- IPCM использует полости на основе поверхности постоянной плотности вокруг молекулы;
- SCIPCM делает то же самое, что и IPCM, только итерационно с изменением поляризации;
метод для учёта влияния электролитов CPCM.

Термодинамическое моделирование

С точки зрения термодинамики любая химическая система стремится к состоянию равновесия, которое соответствует минимуму свободной энергии системы. Термодинамическое моделирование заключается в нахождении выражения для свободной энергии системы при заданных условиях, поиск ее минимума, – точки, которая определяет равновесный состав системы (или другие искомые параметры).

В этой области существуют различные программные продукты, в том числе сервис SlobodovReactor. Подробнее об этом методе можно прочитать в книге «Компьютерное моделирование химических взаимодействий».