Gaussian-1 и Gaussian-2: обобщенные методы
для предсказания энергии молекул

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

Gaussian-1

Pople с сотрудниками предложили такой общий метод, пригодный практически для любых молекулярных систем, позволяющий предсказывать экспериментальные результаты с заданной точностью. Для таких процессов как диссоциация, ионизация и простые реакции различие с экспериментом не превышает 2 ккал/моль. Этот метод авторы назвали Gaussian-1 или G1.

Разрабатывая метод G1 авторы попытались преодолеть недостатки простых методов. Например, известно, что однодетерминантный метод Хартри-Фока плохо воспроизводит энергии диссоциации связей, если не проводится учет энергии корреляции электронов. В G1 электронная корреляция учитывается с помощью процедуры Меллера-Плессета вплоть до 4-го порядка. Однако точность достигаемая с помощью MP обычно недостаточна, поэтому дальнейшая корреляция учитывается методом квадратичного конфигурационного взаимодействия. Далее учитывается вклад d- и f-поляризационных функций. Кроме того, производится поправка на остаточные неточности, возникающие вследствие присущего HF методу неучета отталкивания между спаренными электронами. Для этого использована концепция изогиричности (изоспиновости), т.е. делается поправка к общей энергии с учетом систем H и H2, для которых существуют наиболее точные решения.

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

Итак, назначение метода состоит в получении

  1. Равновесной геометрии
  2. Полной электронной энергии
  3. Набора гармонических частот соответствующих локальному минимуму на потенциальной поверхности.

Полная электронная энергия Ee соответствует энергии системы при фиксированных положениях ядер, соответствующих минимуму. Далее, используя гармонические частоты можно рассчитать коррекцию на нулевую точку DE(ZPE) и, следовательно, энергию низшего колебательного состояния молекулы:

E0 = Ee + DE(ZPE),

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

Теперь можно описать последовательность шагов для определения вышеназванных величин.

Шаг 1. Начальная геометрия рассчитывается однодетерминантным методом Хартри-Фока с использованием базиса 6-31G* (HF/6-31G*). Применяется ограниченный HF для синглетов (RHF) и неограниченный для состояний с другими мультиплетностями (UHF). В последующих шагах применяется тот же метод, что и в первом (RHF или UHF).

Шаг 2. Равновесная геометрия рассчитывается с использованием метода Меллера-Плессета второго порядка с базисом 6-31G*. Используется MP2 с полной электронной корреляцией (включает электроны с внутренних оболочек). Полученная таким образом геометрия обозначается как MP2(FU)/ 6-31G*.

Шаг 3. Структура, полученная на шаге 2, используется в ряде расчетов с фиксированной геометрией с применением более высоких уровней электронной теории. Все они используют электронную корреляцию только для валентных электронов (FC приближение).

Первый ряд высокоуровневых вычислений проводится полным методом Меллера-Плессета четвертого порядка (MP4SDTQ) с базисом 6-311G**. Полное обозначение этого метода

MP4SDTQ(FC)/6-311G**//MP2(FU)/6-31G*

.Далее FC и обозначение геометрии будет опускаться.

MP4SDTQ/6-311G** энергии еще недостаточно точны, чтобы дать удовлетворительное значение Ee, поэтому они подвергаются ряду поправок. Для получения этих поправок необходимо проведение довольно значительного объема расчетов. Поэтому поправки рассчитываются отдельно, а затем суммируются, при этом предполагается, что они аддитивны.

Шаг 4. Первая поправка включает учет диффузных sp базисных функций, которые, как известно, важны для ионов и молекул с неподеленными парами электронов:

DE(+) = E(MP4/6-311+G**) – E(MP4/6-311G**)

Шаг 5. Вторая поправка включает учет высших поляризационных функций на неводородных атомах:

DE(2df) = E(MP4/6-311G**(2df)) – E(MP4/6-311G**)

Шаг 6. Третья поправка учитывает недостаток метода Меллера-Плессета, усеченного до четвертого порядка. Известно, что иногда MP медленно выходит на предел при расчете молекул с напряженными связями. Однако сходимость бывает плохой и при равновесных связях, и ограничение четвертым порядком может приводить к значительным ошибкам. В частности, это применимо для ненасыщеных радикалов и молекул с кратными связями.

Для преодоления вышеназванных трудностей, далее применяется пост-MP4 поправка на остаточные корреляционные эффекты. Она основана на теории квадратичного конфигурационного взаимодействия QCISD(T), которая хорошо воспроизводит полное конфигурационное взаимодействие (FCI), в частности около равновесной геометрии. Метод QCISD(T) эквивалентен MP4 в четвертом порядке, а также учитывает поправки пятого и высших порядков. Таким образом, поправка на квадратичное конфигурационное взаимодействие представляет собой:

DE(QCI) = E(QCISD(T)/6-311G**) – E(MP4/6-311G**)

Складывая три полученных поправки и добавляя их к E(MP4/6-311G**) мы получаем комбинированную энергию:

E(combined) = E(MP4/6-311G**) + DE(+) + DE(2df) + DE(QCI),

которая представляет собой приближение к полному расчету на уровне QCISD(T)/6-311+G**(2df). Таким образом, один большой расчет мы заменили на четыре значительно меньших.

Шаг 7. Энергии, полученные на последнем шаге все еще неполны, в основном вследствие неполноценности базиса. Главная проблема состоит в том, что корреляционная энергия двух спаренных электронов очень медленно выходит на предел с расширением базиса. Этот эффект исправляется не проведением дополнительных квантовохимических расчетов, а наложением общей высокоуровневой поправки DE(HLC), которая приводит к верным значениям Ee для атома и молекулы водорода. Добавляется -3.07 миллихартри на каждый спаренный валентный электрон и -0.19 миллихартри на каждый неспаренный валентный электрон. Если мы рассматриваем спиновое состояние с максимальным числом a-электронов, то полная коррекция (в миллихартри) составит:

DE(HLC) = – 0.19na – 5.95nb,

где na и nb представляют собой число a и b валентных электронов соответственно.

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

Шаг 7 дает следующее значение электронной энергии:

Ee = E(combined) + DE(HLC)

Шаг 8. Последний шаг включает расчет гармонических частот и добавление коррекции на нулевую температуру DE(ZPE) с получением полной энергии E0. Частоты рассчитываются с использованием уровня 6-31G* (используя Хартри-Фоковскую структуру) и умножаются на фактор 0.8929 для учета неадекватности воспроизведения частот на этом уровне теории.

Gaussian-2

G1 плохо воспроизводит энергии диссоциации ионных молекул (LiF), триплетные состояния молекул (O2, S2), разделение энергии в триплетах и синглетах (CH2, SiH2), некоторые гидриды (NH3, N2H4), гипервалентные молекулы (SO2, ClO2). В процедуре G2 устранен ряд недостатков метода G1, что позволило существенно улучшить получаемые результаты. Введены следующие новые возможности:

  1. Исключено предположение об аддитивности диффузных sp (+) и 2df расширений базиса, которое использовалось в G1. Это привело к существенному улучшению результатов для молекул с ионным характером и анионов.
  2. Добавлены третья d функция для неводородных атомов и вторая p функция для атомов водорода. Третья d функция особенно важна для гипервалентных молекул, содержащих атомы из третьего периода, такие как ClO2 и SO2. Вторая p функция существенно улучшает энергии атомизации для ряда водородсодержащих молекул.
  3. Высокоуровневая поправка (HLC), которая определялась в G1 как ошибка в рассчитанной энергии молекулы и атома водорода, теперь выведена из экспериментальных энергий атомизации для 55 молекул.

Для определения энергии G2 необходимо вычислить поправку D к G1 энергии. Она состоит из двух частей рассчитываемых на MP2 уровне. Первая поправка D1 исправляет неаддитивность, вызванную разделением расширения базиса для диффузных sp функций (+) и высших поляризационных функций (2df):

D1 = D(+2df) – D(+) – D(2df),

где D(+2df) = E(MP2/6-311+G(2df,p)) – E(MP2/6-311G(d,p)),

D(+) = E(MP2/6-311+G(d,p)) – E(MP2/6-311G(d,p)) и

D(2df) = E(MP2/6-311G(2df,p)) – E(MP2/6-311G(d,p)).

Вторая часть D2 представляет поправку на добавление третьей d функции на неводородные атомы и второй p функции на атомы водорода:

D2 = E(MP2/6-311+G(3df,2p)) – E(MP2/6-311+G(2df,p)).

Суммарная поправка D тогда будет равна

D = D1 + D2 = E(MP2/6-311+G(3df,2p))

E(MP2/6-311G(2df,p))

E(MP2/6-311+G(d,p))

+ E(MP2/6-311G(d,p))

Таким образом, G2 метод требует проведения только одного дополнительного расчета E(MP2/6-311+G(3df,2p)), а остальные MP2 значения уже получены G1 методом как часть MP4 расчетов.

G2 энергия тогда определяется следующим выражением:

E0(G2) = E0(G1) + D + 1.14 npair ,

где npair – это число валентных пар. Последнее слагаемое, понижающее высокоуровневую поправку HLC в G1 на 1.14 миллихартри, получено из наиболее надежных экспериментальных значений энергии атомизации для 55 молекул.

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

Gaussian-2M

Mebel с сотрудниками предложили серию приближений к G2, названных G2M. Расчет колебательных частот для определения нулевой колебательной энергии на уровне HF и оптимизация геометрии с учетом корреляции электронов MP2 заменены на расчет по теории функционала плотности (DFT) с использованием функционала B3LYP. QCISD(T) расчеты заменены на ограниченный по спину CCSD(T), что еще улучшило точность метода. Это приближение выбрали также Bauschlicher и Partridge, получившие сходную процедуру. Ожидалось, что G2M(RCC) будет хорошо работать в системах со значительным спиновым загрязнением волновой функции нулевого порядка, что обычно для переходных состояний.