Теория возмущений Мёллера–Плессета (MP)

Теория возмущений Мёллера–Плессета (MP) – это метод учёта динамической корреляции, являющийся альтернативой методу конфигурационного взаимодействия.

Примечание

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

В квантовой химии применяется теория возмущения Рэлея–Шрёдингера (RS) – разновидность теории возмущений, основанная на разложении энергии и волновой функции в ряды Тейлора.

Также данный метод часто называется теорией возмущения Мёллера–Плессета (MP) в честь первых учёных, применивших теорию возмущения Рэлея–Шрёдингера к многоэлектронным системам.

Ещё одно название данного метода – теория возмущения многих тел (MBPT) – наиболее часто применяется в контексте бесконечных периодических систем (например, кристаллов).

В основе метода лежит следующая идея: для поиска выражения точной волновой функции \(\Psi\) необходимо получить набор функций, ортогональных по отношению друг к другу и по отношению к волновой функции \(\Psi_0\), полученной методом HF. После этого полученный набор функций вместе с функцией \(\Psi_0\) можно использовать в качестве базисного набора для поиска точной волновой функции \(\Psi\). В качестве такого набора используется набор возбуждённых конфигураций \(\Psi_i\), полученных из основной конфигурации \(\Psi_0\):

\[\begin{split}\Psi= c_{0}\Psi_{0}+ \sum_{\substack{a\\r}}c^{r}_{a}\Psi^{r}_{a}+ \sum_{\substack{a<b\\r<s}}c^{rs}_{ab}\Psi^{rs}_{ab}+ \sum_{\substack{a<b<c\\r<s<t}}c^{rst}_{abc}\Psi^{rst}_{abc}+ \sum_{\substack{a<b<c<d\\r<s<t<u}}c^{rstu}_{abcd}\Psi^{rstu}_{abcd}+ \ldots\end{split}\]

Теория возмущений предлагает альтернативный подход: волновые функции \(\Psi_i\) представляются в виде суммы нулевого приближения \(\Psi^{(0)}_i\) и поправок \(\Psi^{(n\ne0)}_i\):

\[\Psi_i=\Psi^{(0)}_i+\Psi^{(1)}_i+\Psi^{(2)}_i+\Psi^{(3)}_i+\ldots\]

Но как же находятся данные поправки?

Рассмотрим теорию возмущений несколько подробнее:

Важно

Прежде чем перейти к подробному рассмотрению необходимо определиться с используемыми понятиями:

Что такое электронное состояние системы?

В стационарных условиях, то есть, если система не зависит от времени, электронное состояние системы – это полный набор характеризующих её параметров, соответствующих определённому распределению электронов по спин-орбиталям. Если же условия нестационарны, то к данным параметрам добавляются ещё параметры, характеризующие её поведение во времени.

В литературе, особенно часто в посвящённой многоконфигурационным методам, волновую функцию, описывающую электронное состояние системы именуют, просто, состоянием.

Конфигурацией же именуется волновая функция, построенная на основе параметров системы, соответствующих определённому распределению электронов по спин-орбиталям.

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

Поэтому в рамках одноконфигурационных методов, например, метода Хартри-Фока, понятия состояние и конфигурация эквивалентны, но для многоконфигурационных методов это уже, как правило, не так.

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

Один из немногих многоконфигурационных методов, где данные понятия эквивалентны – это метод CIS, так как референтная конфигурация в нём не взаимодействует с остальными, из-за чего она, фактически, равна референтного состоянию.

В основе метода теории возмущений лежит предположение, что точный гамильтониан системы \(\hat{H}\) отличается от уже изветного нулевого приближения \(\hat{H}_0\) на достаточно «малый» оператор возмущения \(\hat{V}\). Таким образом:

\[\hat{H}=\hat{H}_0+\hat{V}\]

Чем же являются нулевое приближение к гамильтониану, оператор возмущения и функции \(\Psi^{(0)}_i\)? В теории возмущения Мёллера-Плессета в качестве нулевого приближения к гамильтониану берётся электронный гамильтониан в методе Хартри-Фока:

\[\hat{H}_0=\sum_{i}\hat{h}(q_i)+\sum_{i}\hat{v}(q_i)\]

Примечание

Так как взаимное отталкивание ядер в методе Хартри-Фока описывается точно (ведь, согласно приближению Борна-Оппенгеймера, оно не зависит от координат электронов), то в теории возмущений Мёллера-Плессета рассматривается только электронный гамильтониан.

Оператор возмущения же выбирается следующим образом:

\[\hat{V}=\hat{H}{(ne)}-\sum_i\hat{v}(q_i)\]

То есть, это оператор, собственным значением которого является ошибка в энергии, возникающая вследствие замены многоэлектронного оператора \(\hat{H}{(ne)}\), описывающего взаимодействие электронов, на сумму эффективных потенциалов \(\sum_i\hat{v}(q_i)\) электронов в усреднённом поле, создаваемом другими электронами.

Иными словами – это оператор ошибки, вносимой методом Хартри-Фока при поиске аппроксимации \(\Psi^{(0)}_0\) к волновой функции \(\Psi\).

В качестве функции \(\Psi^{(0)}_0\) в теории возмущения Мёллера-Плессета выбирается конфигурация \(\Psi_0\), полученная методом Хартри-Фока. А в качестве остальных (\(\Psi^{(0)}_{i\ne0}\)) – возбуждённые конфигурации \(\Psi_i\).

Для упрощения следующих математических преобразований введём безразмерный множитель \(\lambda=1\).

Домножим в выражении гамиольтониана \(\hat{H}\) оператор \(\hat{V}\) на \(\lambda\):

\[\hat{H}=\hat{H}_0+\lambda\hat{V}\]

Домножив функции \(\Psi^{(n)}_i\) на множители \(\lambda^n\), получает выражения волновой функции \(\Psi_i\) в виде ряда Тейлора по множителю \(\lambda\):

\[\Psi_i=\Psi^{(0)}_i+\lambda\Psi^{(1)}_i+\lambda^2\Psi^{(2)}_i+\lambda^3\Psi^{(3)}_i+\ldots\]

Введём аналогичный ряд для энергии \(E_i\):

\[E_i=E^{(0)}_i+\lambda{E}^{(1)}_i+\lambda^2{E}^{(2)}_i+\lambda^3{E}^{(3)}_i+\ldots\]

где \(E^{(0)}_i\) – нулевое приближение к \(E_i\), а \(E^{(1)}_i\), \(E^{(2)}_i\), \(E^{(3)}_i\), \(\ldots\) – поправки к нему.

Для поправок \(\Psi^{(1)}_i\), \(\Psi^{(2)}_i\), \(\Psi^{(3)}_i\), \(\ldots\) ставится условие ортогональнальности нулевому приближению \(\Psi^{(0)}_i\):

\[ \begin{align}\begin{aligned}\int\Psi^{(0)*}_i\Psi^{(1)}_idx=0\\\int\Psi^{(0)*}_i\Psi^{(2)}_idx=0\\\int\Psi^{(0)*}_i\Psi^{(3)}_idx=0\\...................\end{aligned}\end{align} \]

Примечание

Так как:

\[\int\Psi^{(0)*}_i\Psi^{(0)}_idx=1\]

То вышеуказанные условия равносильны:

\[\int\Psi^{(0)*}_i\Psi_idx=1\]

Теперь возьмём уравнение Шрёдингера:

\[\hat{H}\Psi_i=E_i\Psi_i\]

И подставим в него записанные выше выражения для гамильтониана, энергии и волновой функции:

\[\begin{split}(\hat{H}_0+\lambda\hat{V})(\Psi^{(0)}_i+\lambda\Psi^{(1)}_i&+\lambda^2\Psi^{(2)}_i+\lambda^3\Psi^{(3)}_i+\ldots)= \\ =(E^{(0)}_i+\lambda{E}^{(1)}_i+\lambda^2{E}^{(2)}_i+\lambda^3{E}^{(3)}_i+\ldots&)(\Psi^{(0)}_i+\lambda\Psi^{(1)}_i +\lambda^2\Psi^{(2)}_i+\lambda^3\Psi^{(3)}_i+\ldots)\end{split}\]

Раскроем скобки и сгруппируем слагаемые по степеням \(\lambda\):

\[\begin{split}\begin{split} \hat{H}_0\Psi^{(0)}_i+\lambda(\hat{H}_0\Psi^{(1)}_i+\hat{V}\Psi^{(0)}_i)+\lambda^2(\hat{H}_0&\Psi^{(2)}_i+\hat{V}\Psi^{(1)}_i) +\lambda^3(\hat{H}_0\Psi^{(3)}_i+\hat{V}\Psi^{(2)}_i)+\ldots= \\ =E^{(0)}_i\Psi^{(0)}_i+\lambda(E^{(0)}_i\Psi^{(1)}_i+E^{(1)}_i\Psi^{(0)}_i&) +\lambda^2(E^{(0)}_i\Psi^{(2)}_i+E^{(1)}_i\Psi^{(1)}_i+E^{(2)}_i\Psi^{(0)}_i)+ \\ +\lambda^3(E^{(0)}_i\Psi^{(3)}_i+E^{(1)}_i\Psi^{(2)}_i&+E^{(2)}_i\Psi^{(1)}_i+E^{(3)}_i\Psi^{(0)}_i)+\ldots \end{split}\end{split}\]

Приравняв слагаемые при одинаковых степенях \(\lambda\) в левой и правой частях, получим систему уравнений:

\[\begin{split}\left\{ \begin{aligned} &\hat{H}_0\Psi^{(0)}_i\phantom{\hat{V}+\Psi^{(0)}_i}=E^{(0)}_i\Psi^{(0)}_i \\ &\hat{H}_0\Psi^{(1)}_i+\hat{V}\Psi^{(0)}_i=E^{(0)}_i\Psi^{(1)}_i+E^{(1)}_i\Psi^{(0)}_i \\ &\hat{H}_0\Psi^{(2)}_i+\hat{V}\Psi^{(1)}_i=E^{(0)}_i\Psi^{(2)}_i+E^{(1)}_i\Psi^{(1)}_i+E^{(2)}_i\Psi^{(0)}_i \\ &\hat{H}_0\Psi^{(3)}_i+\hat{V}\Psi^{(2)}_i=E^{(0)}_i\Psi^{(3)}_i+E^{(1)}_i\Psi^{(2)}_i+E^{(2)}_i\Psi^{(1)}_i+E^{(3)}_i\Psi^{(0)}_i \\ &........................................................... \end{aligned} \right.\end{split}\]

Примечание

Как видно из первого выражения, все конфигурации \(\Psi^{(0)}_i\) являются собственными функциями оператора \(\hat{H}_0\).

Каждое из полученных уравнений домножим слева на \(\Psi^{(0)*}_i\) и проинтегрируем по \(x\). Приняв во внимание упомянутые выше условия ортогональности, получим:

\[\begin{split}\left\{ \begin{aligned} &E^{(0)}_i=\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(0)}_idx \\ &E^{(1)}_i=\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(1)}_idx+\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_idx \\ &E^{(2)}_i=\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(2)}_idx+\int\Psi^{(0)*}_i\hat{V}\Psi^{(1)}_idx \\ &E^{(3)}_i=\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(3)}_idx+\int\Psi^{(0)*}_i\hat{V}\Psi^{(2)}_idx \\ &.......................................... \end{aligned} \right.\end{split}\]

Примечание

Докажем, что слагаемые:

\[\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(n)}_idx\]

Равны нулю при \(n\ne0\).

Запишем уравнение Шрёдингера:

\[\hat{H}_0\Psi^{(0)}_i=E^{(0)}_i\Psi^{(0)}_i\]

Домножим его на \(\Psi^{(n)*}_i\) и проинтегрируем по \(x\):

\[\int\Psi^{(n)*}_i\hat{H}_0\Psi^{(0)}_idx=\int\Psi^{(n)*}_iE^{(0)}_i\Psi^{(0)}_idx\]

Вынесем \(E^{(0)}_i\) из под интеграла:

\[\int\Psi^{(n)*}_i\hat{H}_0\Psi^{(0)}_idx=E^{(0)}_i\int\Psi^{(n)*}_i\Psi^{(0)}_idx\]

Так как оператор \(\hat{H}_0\) – эрмитов, то справедливы следующие преобразования:

\[ \begin{align}\begin{aligned}\left(\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(n)}_idx\right)^*&=E^{(0)}_i\left(\int\Psi^{(0)*}_i\Psi^{(n)}_idx\right)^*\\\left(\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(n)}_idx\right)^*&=\left(E^{(0)}_i\int\Psi^{(0)*}_i\Psi^{(n)}_idx\right)^*\\\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(n)}_idx&=E^{(0)}_i\int\Psi^{(0)*}_i\Psi^{(n)}_idx\end{aligned}\end{align} \]

Как было отмечено выше, по условию ортогональности интеграл в правой части при \(n\ne0\) равен нулю, а значит, интеграл в левой части также равен нулю.

Так как все слагаемые вида:

\[\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(n)}_idx\]

Равны нулю при \(n\ne0\), то система уравнений упрощается и мы получаем выражения для всех \(E^{(n)}_i\):

\[\begin{split}\left\{ \begin{aligned} &E^{(0)}_i=\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(0)}_idx \\ &E^{(1)}_i=\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_idx \\ &E^{(2)}_i=\int\Psi^{(0)*}_i\hat{V}\Psi^{(1)}_idx \\ &E^{(3)}_i=\int\Psi^{(0)*}_i\hat{V}\Psi^{(2)}_idx \\ &............................ \end{aligned} \right.\end{split}\]

Но вернёмся к исходной системе уравнений и возьмём из неё второе уравнение:

\[\hat{H}_0\Psi^{(1)}_i+\hat{V}\Psi^{(0)}_i=E^{(0)}_i\Psi^{(1)}_i+E^{(1)}_i\Psi^{(0)}_i\]

Перенесём слагаемые, содержащие \(\Psi^{(1)}_i\), в левую часть, а остальные – в правую, после чего вынесем общие множители за скобку и, для удобства, умножим обе части на \(-1\):

\[ \begin{align}\begin{aligned}\hat{H}_0\Psi^{(1)}_i-E^{(0)}_i\Psi^{(1)}_i&=E^{(1)}_i\Psi^{(0)}_i-\hat{V}\Psi^{(0)}_i\\(\hat{H}_0-E^{(0)}_i)\Psi^{(1)}_i&=(E^{(1)}_i-\hat{V})\Psi^{(0)}_i\\(E^{(0)}_i-\hat{H}_0)\Psi^{(1)}_i&=(\hat{V}-E^{(1)}_i)\Psi^{(0)}_i\end{aligned}\end{align} \]

Теперь представим \(\Psi^{(1)}_i\) в виде разложения по базису из конфигураций \(\Psi^{(0)}_j\):

\[(E^{(0)}_i-\hat{H}_0)\sum_{j}c^{(1)}_j\Psi^{(0)}_j=(\hat{V}-E^{(1)}_i)\Psi^{(0)}_i\]

Так как по условию функции \(\Psi^{(0)}_i\) и \(\Psi^{(1)}_i\) ортогональны, то коэффициент в разложении \(c^{(1)}_i\) будет равен нулю. Поэтому соответствующее ему слагаемое можно исключить из суммы. Сумму с исключённым слагаемым пометим штрихом (\('\)):

\[(E^{(0)}_i-\hat{H}_0)\sum_{j}'c^{(1)}_j\Psi^{(0)}_j=(\hat{V}-E^{(1)}_i)\Psi^{(0)}_i\]

Теперь домножим обе части на \(\Psi^{(0)*}_k\) и проинтегрируем по \(x\):

\[\int\Psi^{(0)*}_k(E^{(0)}_i-\hat{H}_0)\sum_{j}'c^{(1)}_j\Psi^{(0)}_jdx=\int\Psi^{(0)*}_k(\hat{V}-E^{(1)}_i)\Psi^{(0)}_idx\]

Примечание

Так как конфигурации \(\Psi^{(0)}_j\) – собственные функции оператора \(\hat{H}_0\), то они ортонормированы, то есть:

\[ \begin{align}\begin{aligned}\int\Psi^{(0)*}_j\Psi^{(0)}_jdx=1\\\int\Psi^{(0)*}_k\Psi^{(0)}_jdx=0\end{aligned}\end{align} \]

При \(j\ne{k}\).

Кроме того:

\[ \begin{align}\begin{aligned}\int\Psi^{(0)*}_j\hat{H}_0\Psi^{(0)}_jdx&=E^{(0)}_j\\\int\Psi^{(0)*}_k\hat{H}_0\Psi^{(0)}_jdx&=0\end{aligned}\end{align} \]

При \(j\ne{k}\).

Учитывая ортонормированность функций \(\Psi^{(0)}_j\), при раскрытии скобок получим:

\[(E^{(0)}_i-E^{(0)}_k)c^{(1)}_k=\int\Psi^{(0)*}_k\hat{V}\Psi^{(0)}_idx\]

Отсюда, при замене \(k\) на \(j\), получаем выражение для коэффициентов \(c^{(1)}_j\):

\[c^{(1)}_j=\frac{\int\Psi^{(0)*}_j\hat{V}\Psi^{(0)}_idx}{E^{(0)}_i-E^{(0)}_j}\]

Таким образом, волновая функция \(\Psi^{(1)}_i\) имеет вид:

\[\Psi^{(1)}_i=\sum_{j}'\frac{\int\Psi^{(0)*}_j\hat{V}\Psi^{(0)}_idx}{E^{(0)}_i-E^{(0)}_j}\Psi^{(0)}_j\]

Теперь возьмём выражение для энергии \(E^{(2)}_i=\int\Psi^{(0)*}_i\hat{V}\Psi^{(1)}_idx\) и подставим в него выражение для функции \(\Psi^{(1)}_i\):

\[E^{(2)}_i=\int\Psi^{(0)*}_i\hat{V}\sum_{j}'\frac{\int\Psi^{(0)*}_j\hat{V}\Psi^{(0)}_idx}{E^{(0)}_i-E^{(0)}_j}\Psi^{(0)}_jdx\]

Вынесем знак суммы за интеграл:

\[E^{(2)}_i=\sum_{j}'\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_jdx\frac{\int\Psi^{(0)*}_j\hat{V}\Psi^{(0)}_idx}{E^{(0)}_i-E^{(0)}_j}\]

Откуда получаем выражение для энергии второго порядка (\(E^{(2)}_i\)):

\[E^{(2)}_i=\sum_{j}'\frac{\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_jdx\int\Psi^{(0)*}_j\hat{V}\Psi^{(0)}_idx}{E^{(0)}_i-E^{(0)}_j} =\sum_{j}'\frac{\left|\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_jdx\right|^2}{E^{(0)}_i-E^{(0)}_j}\]

Аналогичным образом можно также получить выражения для энергий третьего (\(E^{(3)}_i\)), четвёртого (\(E^{(4)}_i\)) и т. д. порядков.

Если принять для краткости следующие обозначения:

\[ \begin{align}\begin{aligned}\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(0)}_jdx=\langle{i}|{\hat{H}_0}|{j}\rangle\\\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_jdx=\langle{i}|{\hat{V}}|{j}\rangle\end{aligned}\end{align} \]

То выражения для энергий до четвёртого порядка примут следующий вид:

\[\begin{split}E^{(0)}_i&=\langle{i}|{\hat{H}_0}|{i}\rangle \\ E^{(1)}_i&=\langle{i}|{\hat{V}}|{i}\rangle \\ E^{(2)}_i&=\sum_{j}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{i}\rangle} {E^{(0)}_i-E^{(0)}_j} \\ E^{(3)}_i&=\sum_{jk}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)\left(E^{(0)}_i-E^{(0)}_k\right)}- E^{(1)}_i\sum_{j}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)^2} \\ E^{(4)}_i&=\sum_{jkl}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{l}\rangle\langle{l}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)\left(E^{(0)}_i-E^{(0)}_k\right)\left(E^{(0)}_i-E^{(0)}_l\right)}- \\ &-E^{(1)}_i\sum_{jk}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)^2\left(E^{(0)}_i-E^{(0)}_k\right)}- E^{(1)}_i\sum_{jk}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)\left(E^{(0)}_i-E^{(0)}_k\right)^2}+ \\ &+\left(E^{(1)}_i\right)^2\sum_{j}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)^3}- E^{(2)}_i\sum_{j}'\frac{\langle{i}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{i}\rangle} {\left(E^{(0)}_i-E^{(0)}_j\right)^2} \\\end{split}\]

Как видно из выражений энергий \(E^{(n)}_i\) выше, сложность вычислений с помощью теории возмущений Мёллера-Плессета (MP) при переходе к энергиям более высокого порядка возрастает с огромной скоростью, поэтому на практике расчёт ограничивают до энергии \(E^{(n)}_i\), где число \(n\) называют порядком теории. Отсюда исходят и принятые обозначения методов, например, MP2 – это теория возмущений Мёллера-Плессета второго порядка.

Но чему же обычно равно \(n\)? Для ответа на данный вопрос вернёмся к определению электронного гамильтониана:

\[\hat{H}_0=\sum_{j}\hat{h}(q_j)+\sum_{j}\hat{v}(q_j)\]

\[\hat{V}=\hat{H}{(ne)}-\sum_j\hat{v}(q_j)\]

Из этого следует, что:

\[\hat{H}=\hat{H}_0+\hat{V}=\left(\sum_{j}\hat{h}(q_j)+\sum_{j}\hat{v}(q_j)\right)+\left(\hat{H}{(ne)}-\sum_j\hat{v}(q_j)\right)= \sum_{j}\hat{h}(q_j)+\hat{H}{(ne)}\]

А значит:

\[\begin{split}E^{(0)}_i+E^{(1)}_i&=\int\Psi^{(0)*}_i\hat{H}_0\Psi^{(0)}_idx+\int\Psi^{(0)*}_i\hat{V}\Psi^{(0)}_idx= \\ &=\int\Psi^{(0)*}_i(\hat{H}_0+\hat{V})\Psi^{(0)}_idx= \\ &=\int\Psi^{(0)*}_i\left(\sum_{j}\hat{h}(q_j)+\hat{H}{(ne)}\right)\Psi^{(0)}_idx=E^{\mathbf{HF}}_i\end{split}\]

То есть, сумма энергий нулевого \(E^{(0)}_i\) и первого \(E^{(1)}_i\) порядков – это не что иное, как энергия \(E^{\mathbf{HF}}_i\), полученная методом Хартри-Фока.

Но энергии второго \(E^{(2)}_i\), третьего \(E^{(3)}_i\) и т. д. порядков в методе Хартри-Фока не учитываются, а так как сумма энергий всех порядков, по определению, равна точной энергии, то сумма всех энергий от второго порядка и выше соответствует энергии \(E^{\mathbf{corr}}_i\), обусловленной явлением динамической корреляции:

\[ \begin{align}\begin{aligned}&E_i=\underbrace{E^{(0)}_i+{E}^{(1)}_i}+\underbrace{{E}^{(2)}_i+{E}^{(3)}_i+{E}^{(4)}_i+\ldots}\\&\phantom{)))))))))))}E^{\mathbf{HF}}_i\phantom{)))))))))))))))))}E^{\mathbf{corr}}_i\end{aligned}\end{align} \]

Примечание

Стоит отметить, что представление электронного гамильтониана в виде:

\[\hat{H}_0=\sum_{j}\hat{h}(q_j)+\sum_{j}\hat{v}(q_j)\]

Используется в методе Хартри-Фока только в процессе итерационного нахождения аппроксимации к волновой функции. Итоговое значение энергии рассчитывается с помощью классического электронного гамильтониана:

\[\hat{H}=\sum_{j}\hat{h}(q_j)+\hat{H}{(ne)}\]

Почему именно так? – Дело в том, что классический электронный гамильтониан \(\hat{H}\) не позволяет найти аппроксимацию волновой функции методом самосогласованного поля, но, так как он ближе к точному электронному гамильтониану, чем \(\hat{H}_0\), то энергия, полученная с его помощью ближе к точной энергии.

Итак, в пределе, теория возмущений даёт полный учёт динамической корреляции. Но насколько велик вклад каждой компоненты \(E^{(n)}_i\)?

В первую очередь, нас интересует энергия референтного состояния \(E_0\), поэтому проанализируем вклад компонент энергии именно на её примере:

Рассмотрим выражение:

\[\int\Psi^*_0\hat{H}\Psi^r_adx\]

По теореме Бриллюэна данное выражение равно нулю (так как референтная и однократно возбуждённые конфигурации не взаимодействуют).

С другой стороны:

\[\int\Psi^*_0\hat{H}\Psi^r_adx=\int\Psi^*_0\left(\hat{H}_0+\hat{V}\right)\Psi^r_adx=\int\Psi^*_0\hat{H}_0\Psi^r_adx+\int\Psi^*_0\hat{V}\Psi^r_adx\]

Рассмотрим первое слагаемое:

Электронный гамильтониан \(\hat{H}_0\) можно представить в виде суммы одноэлектронных операторов Фока:

\[\hat{H}_0=\sum_{i}\hat{f}(x_i)\]

Тогда:

\[\int\Psi^*_0\hat{H}_0\Psi^r_adx=\int\Psi^*_0\sum_{i}\hat{f}(x_i)\Psi^r_adx=\sum_{i}\int\Psi^*_0\hat{f}(x_i)\Psi^r_adx\]

Теперь рассмотрим конфигурации: их можно представить в виде определителя Слэтера. Но что из себя представляет определитель Слэтера? – Для конфигурации \(\Psi_0\), по определению, он равен:

\[\Psi_0=\frac{1}{\sqrt{n!}}\sum^{n!}_{\alpha=1}(-1)^{p_{\alpha}}\hat{P}_{\alpha}\psi_1(x_1)\psi_2(x_2)\ldots\psi_n(x_n)\]

где:

• \(\hat{P}_{\alpha}\) – оператор перестановки: данный оператор переставляет спин-орбитали в произведении после него (порядок координат электронов сохраняется). Так как всего орбиталей \(n\), то число перестановок равно \(n!\).

• \(p_{\alpha}\) – число элементарных (попарных) перестановок, эквивалентных действию оператора перестановки \(\hat{P}_{\alpha}\).

При подстановке его в выражение выше множители \(\frac{1}{\sqrt{n!}}\) и \((-1)^{p_{\alpha}}\), так как они не зависят от координат электронов, всегда можно будет вынести за знаки гамильтониана и интеграла. Поэтому, является ли интеграл нулевым или нет, определяет только произведение спин-орбиталей.

Рассмотрим произвольное произведение спин-орбитали из данной суммы:

\[\psi_k(x_1)\psi_l(x_2)\ldots\psi_z(x_{i-1})\psi_a(x_i)\psi_b(x_{i+1})\ldots\psi_m(x_n)\]

Для конфигурации \(\Psi^r_a\) оно имеет аналогичный вид, единственное различие: спин-орбиталь \(\psi_a\) заменена на спин-орбиталь \(\psi_r\):

\[\psi_k(x_1)\psi_l(x_2)\ldots\psi_z(x_{i-1})\psi_r(x_i)\psi_b(x_{i+1})\ldots\psi_m(x_n)\]

При разложении конфигураций в виде определителя Слэтера выражение \(\sum^n_{i=1}\int\Psi^*_0\hat{f}(x_i)\Psi^r_adx\) обратится в сумму из \(n\cdot(n!)^2\) слагаемых, при оценке значений которых возможны три случая:

1.) Одинаковые спин-орбитали в конфигурациях \(\Psi_0\) и \(\Psi^r_a\) заселяют одинаковые электроны. Спин-орбитали \(\psi_a\) и \(\psi_r\) заселяет \(i\)-ый электрон.

2.) Одинаковые спин-орбитали в конфигурациях \(\Psi_0\) и \(\Psi^r_a\) заселяют одинаковые электроны. Спин-орбитали \(\psi_a\) и \(\psi_r\) заселяет \(j\)-ый электрон (\(i{\ne}j\)).

3.) Одинаковые спин-орбитали в конфигурациях \(\Psi_0\) и \(\Psi^r_a\) заселяют разные электроны.

Для сокрацения записи примем обозначение:

\[\psi_q(x_t):=\psi^t_q\]

В первом случае имеем:

\[\int(\psi^1_k\psi^2_l\ldots\psi^{i-1}_z\psi^i_a\psi^{i+1}_b\ldots\psi^n_m)^*\hat{f}(x_i)\psi^1_k\psi^2_l\ldots\psi^{i-1}_z\psi^i_r\psi^{i+1}_b\ldots\psi^n_mdx=\]

\[=\int\psi^{i*}_a\hat{f}(x_i)\psi^i_rdx_i\int(\psi^1_k\psi^2_l\ldots\psi^{i-1}_z\psi^{i+1}_b\ldots\psi^n_m)^*\psi^1_k\psi^2_l\ldots\psi^{i-1}_z\psi^{i+1}_b\ldots\psi^n_mdx'\]

где \(dx'=dx_1dx_2{\ldots}dx_{i-1}dx_{i+1}{\ldots}dx_n\).

Второй множитель из полученного выражения разбивается на \(n-1\) множителей вида:

\[\int\psi^{t*}_q\psi^t_qdx_t\]

Так как спин-орбитали ортонормированы, то все данные множители равны единице. Тогда:

\[\int(\psi^1_k\psi^2_l\ldots\psi^{i-1}_z\psi^i_a\psi^{i+1}_b\ldots\psi^n_m)^*\hat{f}(x_i)\psi^1_k\psi^2_l\ldots\psi^{i-1}_z\psi^i_r\psi^{i+1}_b\ldots\psi^n_mdx= \int\psi^{i*}_a\hat{f}(x_i)\psi^i_rdx_i\]

Теперь рассмотрим второй случай:

\[\int(\psi^1_k\psi^2_l\ldots\psi^{j-1}_z\psi^j_a\psi^{j+1}_b\ldots\psi^n_m)^*\hat{f}(x_i)\psi^1_k\psi^2_l\ldots\psi^{j-1}_z\psi^j_r\psi^{j+1}_b\ldots\psi^n_mdx=\]

\[=\int\psi^{j*}_a\psi^j_rdx_j\int(\psi^1_k\psi^2_l\ldots\psi^{j-1}_z\psi^{j+1}_b\ldots\psi^n_m)^*\hat{f}(x_i)\psi^1_k\psi^2_l\ldots\psi^{j-1}_z\psi^{j+1}_b\ldots\psi^n_mdx'\]

где \(dx'=dx_1dx_2{\ldots}dx_{j-1}dx_{j+1}{\ldots}dx_n\).

Так как спин-орбитали ортонормированы, то первый множитель в выражении выше обращается в ноль. Тогда:

\[\int(\psi^1_k\psi^2_l\ldots\psi^{j-1}_z\psi^j_a\psi^{j+1}_b\ldots\psi^n_m)^*\hat{f}(x_i)\psi^1_k\psi^2_l\ldots\psi^{j-1}_z\psi^j_r\psi^{j+1}_b\ldots\psi^n_mdx=0\]

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

\[\int\psi^{t*}_q\psi^t_sdx_t\]

Который, ввиду ортонормированности спин-орбиталей, равен нулю.

То есть, слагаемое не равно нулю только в первом случае.

Так как в конфигурациях \(\Psi_0\) и \(\Psi^r_a\) одинаковые орбитали заселяют одинаковые электроны, то операторы \(\hat{P}_{\alpha}\) одинаковы, а значит, равны и \(p_{\alpha}\). Итого имеем, что ненулевые слагаемые в выражении определителя будут иметь вид:

\[\int\frac{1}{\sqrt{n!}}(-1)^{p_{\alpha}}\psi^{i*}_a\hat{f}(x_i)\frac{1}{\sqrt{n!}}(-1)^{p_{\alpha}}\psi^i_rdx_i= \frac{1}{n!}\int\psi^{i*}_a\hat{f}(x_i)\psi^i_rdx_i\]

Слагаемых, удовлетворяющих первому случаю при фиксированном значении индекса \(i\), всего \((n-1)!\), тогда:

\[\int\Psi^*_0\hat{H}_0\Psi^r_adx=\sum^{n}_{i=1}\int\Psi^*_0\hat{f}(x_i)\Psi^r_adx=\sum^{n}_{i=1}\frac{(n-1)!}{n!}\int\psi^{i*}_a\hat{f}(x_i)\psi^i_rdx_i= \int\psi^{i*}_a\hat{f}(x_i)\psi^i_rdx_i\]

Итого:

\[\int\Psi^*_0\hat{H}_0\Psi^r_adx=\int\psi^*_a\hat{f}\psi_rdx\]

Важно

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

Так как вывод даже этого правила весьма объёмен, то вывод остальных правил останется за кадром.

На данный момент повествования, помимо вышеприведённого, важны ещё два правила:

1.) Если конфигурации \(\Psi_i\) и \(\Psi_j\) отличаются более чем одной орбиталью, то:

\[\int\Psi^*_i\hat{F}\Psi_jdx=0\]

где \(\hat{F}\) – любой эрмитов оператор, сводящийся к сумме одноэлектронных операторов, в том числе гамильтониан в методе Хартри-Фока \(\hat{H}_0\).

2.) Если конфигурации \(\Psi_i\) и \(\Psi_j\) отличаются более чем двумя орбиталями, то:

\[\int\Psi^*_i\hat{G}\Psi_jdx=0\]

где \(\hat{G}\) – любой эрмитов оператор, сводящийся к сумме одно- и двухэлектронных операторов, в том числе гамильтониан \(\hat{H}\) и оператор возмущения \(\hat{V}\).

Итак, а теперь вспомним, что все спин-орбитали – это собственные функции одноэлектронного оператора Фока, тогда, ввиду ортонормированности спин-орбиталей:

\[\int\Psi^*_0\hat{H}_0\Psi^r_adx=\int\psi^*_a\hat{f}\psi_rdx=\varepsilon_r\int\psi^*_a\psi_rdx=0\]

Теперь вернёмся к выражению:

\[\int\Psi^*_0\hat{H}\Psi^r_adx=\int\Psi^*_0\hat{H}_0\Psi^r_adx+\int\Psi^*_0\hat{V}\Psi^r_adx\]

Выразим из него второе слагаемое правой части и подставим полученные выше значения:

\[\int\Psi^*_0\hat{V}\Psi^r_adx=\int\Psi^*_0\hat{H}\Psi^r_adx-\int\Psi^*_0\hat{H}_0\Psi^r_adx=0-0=0\]

Итак, для оператора возмущения \(\hat{V}\) применимы правила Слэтера-Кондона, а также аналог теоремы Бриллюэна:

\[\int\Psi^*_0\hat{V}\Psi^r_adx=0\]

То есть, отсеяв слагаемые, обращающиеся в нуль по данным правилам, мы можем установить, какие именно конфигурации входят в каждую компоненту \(E^{(n)}_0\).

Введём, для удобства, следующие обозначения:

\[ \begin{align}\begin{aligned}&|{0}\rangle:=\Psi^{(0)}_0\equiv\Psi_0 &&\langle{0}|:=\Psi^{(0)*}_0\equiv\Psi^*_0\\&|{S_i}\rangle:=\Psi^{(0)}_i\sim\Psi^{r}_{a} &&\langle{S_i}|:=\Psi^{(0)*}_i\sim\Psi^{r*}_{a}\\&|{D_i}\rangle:=\Psi^{(0)}_i\sim\Psi^{rs}_{ab} &&\langle{D_i}|:=\Psi^{(0)*}_i\sim\Psi^{rs*}_{ab}\\&|{T_i}\rangle:=\Psi^{(0)}_i\sim\Psi^{rst}_{abc} &&\langle{T_i}|:=\Psi^{(0)*}_i\sim\Psi^{rst*}_{abc}\\&|{Q_i}\rangle:=\Psi^{(0)}_i\sim\Psi^{rstu}_{abcd} &&\langle{Q_i}|:=\Psi^{(0)*}_i\sim\Psi^{rstu*}_{abcd}\\&|{P_i}\rangle:=\Psi^{(0)}_i\sim\Psi^{rstuv}_{abcde} &&\langle{P_i}|:=\Psi^{(0)*}_i\sim\Psi^{rstuv*}_{abcde}\\&|{H_i}\rangle:=\Psi^{(0)}_i\sim\Psi^{rstuvw}_{abcdef} &&\langle{H_i}|:=\Psi^{(0)*}_i\sim\Psi^{rstuvw*}_{abcdef}\end{aligned}\end{align} \]

То есть, если, например, конфигурация \(\Psi^{(0)}_i\) является однократно вырожденной по отношению к \(|{0}\rangle\), то есть \(\Psi^{r}_{a}\), то переобозначим её как \(|{S_i}\rangle\)

Для произвольных конфигураций примем обозначения:

\[ \begin{align}\begin{aligned}&|{i}\rangle:=\Psi^{(0)}_i &&\langle{i}|:=\Psi^{(0)*}_i\\&|{j}\rangle:=\Psi^{(0)}_j &&\langle{j}|:=\Psi^{(0)*}_j\\&|{k}\rangle:=\Psi^{(0)}_k &&\langle{k}|:=\Psi^{(0)*}_k\\&|{l}\rangle:=\Psi^{(0)}_l &&\langle{l}|:=\Psi^{(0)*}_l\\&|{m}\rangle:=\Psi^{(0)}_m &&\langle{m}|:=\Psi^{(0)*}_m\\&|{n}\rangle:=\Psi^{(0)}_n &&\langle{n}|:=\Psi^{(0)*}_n\end{aligned}\end{align} \]

Проанализируем полученные ранее выражения для \(E^{(n)}_i\) (в нашем случае \(i=0\)):

0 и 1) Компоненты \(E^{(0)}_0\) и \(E^{(1)}_0\) содержат, соответственно, выражения \(\langle{0}|{\hat{H}_0}|{0}\rangle\) и \(\langle{0}|{\hat{V}}|{0}\rangle\), которые не равны нулю.

2) Компонента \(E^{(2)}_0\) содержит выражение:

\[\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{0}\rangle\]

где \(j\ne0\).

Согласно правилам Слэтера-Кондона, чтобы данное выражение было не нулевым, \(|{i}\rangle\) не может быть более, чем двукратно возбуждённой по отношению к \(|{0}\rangle\), а по аналогу теоремы Бриллюэна она не должна быть однократно возбуждённой. То есть, \(|{j}\rangle=|{D_j}\rangle\):

\[\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{0}\rangle=\langle{0}|{\hat{V}}|{D_j}\rangle\langle{D_j}|{\hat{V}}|{0}\rangle\]

3) Компонента \(E^{(3)}_0\) содержит выражение:

\[\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{0}\rangle\]

где \(j,k\ne0\).

По аналогии с \(E^{(2)}_0\), чтобы данное выражение не было нулевым, \(|{j}\rangle\) должна равняться \(|{D_j}\rangle\), а \(|{k}\rangle\) – \(|{D_k}\rangle\). При этом \(|{D_j}\rangle\) и \(|{D_k}\rangle\) должны быть таковы, чтобы выражение \(\langle{D_j}|{\hat{V}}|{D_j}\rangle\) не было равно нулю. Значит:

\[\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{0}\rangle= \langle{0}|{\hat{V}}|{D_j}\rangle\langle{D_j}|{\hat{V}}|{D_k}\rangle\langle{D_k}|{\hat{V}}|{0}\rangle\]

Примечание

Также \(E^{(3)}_0\) содержит второе слагаемое, на выражение в котором накладываются те же ограничения, что и для \(E^{(2)}_0\).

4) Компонента \(E^{(4)}_0\) содержит выражение:

\[\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{l}\rangle\langle{l}|{\hat{V}}|{0}\rangle\]

где \(j,k,l\ne0\).

Ограничения на \(|{j}\rangle\) и \(|{l}\rangle\) здесь аналогичны тем, что были в предыдущих случаях: \(|{j}\rangle=|{D_j}\rangle\) и \(|{l}\rangle=|{D_l}\rangle\). запишем промежуточный результат:

\[\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{l}\rangle\langle{l}|{\hat{V}}|{0}\rangle= \langle{0}|{\hat{V}}|{D_j}\rangle\langle{D_j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{D_l}\rangle\langle{D_l}|{\hat{V}}|{0}\rangle\]

где \(k\ne0\).

А теперь рассмотрим конфигурацию \(|{k}\rangle\): согласно условию \(j,k,l\ne0\) и правилам Слэтера-Кондона она может быть однократно (\(|{S_k}\rangle\)), двукратно (\(|{D_k}\rangle\)), трёхкратно (\(|{T_k}\rangle\)) или четырёхкратно (\(|{Q_k}\rangle\)) возбуждённой по отношению к \(|{0}\rangle\), тогда:

\[\begin{split}\langle{0}|{\hat{V}}|{j}\rangle\langle{j}|{\hat{V}}|{k}\rangle\langle{k}|{\hat{V}}|{l}\rangle\langle{l}|{\hat{V}}|{0}\rangle= \langle{0}|{\hat{V}}|{D_j}\rangle\langle{D_j}|{\hat{V}}| \begin{matrix} Q_k \\ T_k \\ D_k \\ S_k \\ \end{matrix} \rangle\langle \begin{matrix} Q_k \\ T_k \\ D_k \\ S_k \\ \end{matrix} |{\hat{V}}|{D_l}\rangle\langle{D_l}|{\hat{V}}|{0}\rangle\end{split}\]

Примечание

Здесь запись в форме столбца используется для обобщённой записи четырёх возможных случаев.

Примечание

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

5, 6, …) Если подобные рассуждения продолжить, то для \(E^{(5)}_0\), \(E^{(6)}_0\), … можно получить выражения:

\[\begin{split}\langle{0}|{\hat{V}}|{D_j}\rangle\langle{D_j}|{\hat{V}}| \begin{matrix} Q_k \\ T_k \\ D_k \\ S_k \\ \end{matrix} \rangle\langle \begin{matrix} Q_k \\ T_k \\ D_k \\ S_k \\ \end{matrix} |{\hat{V}}| \begin{matrix} Q_l \\ T_l \\ D_l \\ S_l \\ \end{matrix} \rangle\langle \begin{matrix} Q_l \\ T_l \\ D_l \\ S_l \\ \end{matrix} |{\hat{V}}|{D_m}\rangle\langle{D_m}|{\hat{V}}|{0}\rangle\end{split}\]

\[\begin{split}\langle{0}|{\hat{V}}|{D_j}\rangle\langle{D_j}|{\hat{V}}| \begin{matrix} Q_k \\ T_k \\ D_k \\ S_k \\ \end{matrix} \rangle\langle \begin{matrix} Q_k \\ T_k \\ D_k \\ S_k \\ \end{matrix} |{\hat{V}}| \begin{matrix} H_l \\ P_l \\ Q_l \\ T_l \\ D_l \\ S_l \\ \end{matrix} \rangle\langle \begin{matrix} H_l \\ P_l \\ Q_l \\ T_l \\ D_l \\ S_l \\ \end{matrix} |{\hat{V}}| \begin{matrix} Q_m \\ T_m \\ D_m \\ S_m \\ \end{matrix} \rangle\langle \begin{matrix} Q_m \\ T_m \\ D_m \\ S_m \\ \end{matrix} |{\hat{V}}|{D_n}\rangle\langle{D_n}|{\hat{V}}|{0}\rangle\end{split}\]

\[.....................................................................\]

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

Важно

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

Например, метод CISD(T) представляет собой метод CISDT, в котором вместо полноценного учёта трёхкратно возбуждённых конфигураций используется часть вклада \(E^{(4)}_0\), содержащая трёхкратно возбуждённые конфигурации. Разумеется, данные вклады неэквивалентны, но они близки, а расчёт через теорию возмущений требует на порядок меньше ресурсов, что делает замену очень выгодной на практике.

Итак, как мы видим, при увелечении порядка теории \(n\) вклад конфигураций с разным порядком возбуждения, в отличие от метода CI, меняется не плавно, а ступенчато, то есть конфигурации с «новыми» порядками возбуждения появляются только при переходе от теорий нечётного порядка к теориям чётного порядка.

Например, при переходе от MP3 (учитываются компоненты до \(E^{(3)}_0\)) к MP4 (до \(E^{(4)}_0\)) добавляется учёт конфигураций \(|{S_i}\rangle\), \(|{T_i}\rangle\) и \(|{Q_i}\rangle\), а при переходе от MP4 к MP5 новые конфигурации не добавляются; из принципиально нового появляется только множитель вида:

\[\begin{split}\langle \begin{matrix} Q_i \\ T_i \\ D_i \\ S_i \\ \end{matrix} |{\hat{V}}| \begin{matrix} Q_j \\ T_j \\ D_j \\ S_j \\ \end{matrix} \rangle\end{split}\]

который «связывает» конфигурации \(|{S_i}\rangle\), \(|{D_i}\rangle\), \(|{T_i}\rangle\) и \(|{Q_i}\rangle\).

Помимо этого, кроме «главной суммы слагаемых» – той, что содержит конфигурации с наибольшими порядками возбуждения в выражениях компонент \(E^{(n)}_0\) присутствуют и другие, причём часть из них входят в итоговые выражения со знаком «минус».

Учитывая, что вклад конфигураций с ростом порядка возбуждения относительно \(|{0}\rangle\) падает, то вклад в энергию \(E_0\) компонент \(E^{(n)}_0\) при увеличении \(n\) меняется неравномерно. Более того: величина \(E^{(n)}_0\), может принимать как отрицательные, так и положительные значения в зависимости от \(n\).

Как следствие, метод MP, в отличие от метода CI, является невариационным методом, иными словами, переход к теории более высокого порядка не всегда приводит к повышению точности результата. В частности, даже возможны, на первый взгляд, парадоксальные результаты, когда энергии, полученные с помощью теорий возмущений некоторых порядков оказываются ниже, чем \(E_0\).

Однако, теория возмущений обладает двумя неоспоримыми преимуществами перед методом CI:

Во-первых, теория возмущений Мёллера-Плессета любого порядка – это неитерационный метод, то есть в нём отсутствует итерационный поиск решения (за рамками метода HF, с помощью которого получается исходная волновая функция \(\Psi^{(0)}_0\)), это делает рассчёт быстрее и позволяет избежать ошибки итерационного метода.

Во-вторых, теория возмущений Мёллера-Плессета – это размерно-согласованный метод, что делает этот метод (теоретически) применимым к системам любого размера.

Итого, на практике, в связи с вышеописанными особенностями, применяется только несколько вариантов теории возмущений:

• MP2 – теория возмущений Мёллера-Плессета второго порядка: данный метод является самым дешёвым в вычислительном плане среди post-HF методов (но на порядок дороже HF). Он учитывает большую часть динамической корреляции (обычно, порядка 80 % вносимой ей энергии) и всегда даёт более точные результаты, чем HF. Часто применяется как компонент в многореферентных методах, а также иногда – в теории функционала плотности (DH-DFT)

• MP4 – значительно более затратный (дороже итерационного шага методов CISD и CCSD), но более точный: учитывает, в среднем, 95 % \(E^{\mathbf{corr}}\). Чаще всего, используется только часть его компоненты \(E^{(4)}_0\) в составе методов CISD(T), CCSD(T), BD(T), а также их модификаций.

• MP5 – ещё более затратный (стоимость на уровне итерационного шага методов CISDT и CCSDT) и примерно такой же по точности, как и MP4. Используется, как правило, только части его компонент \(E^{(4)}_0\) и \(E^{(5)}_0\) в составе таких относительно редко используемых методов, как CCSD(TQ).

В GAMESS (US) в чистом виде представлена только теория возмущений Мёллера-Плессета второго порядка (MP2).

Для того чтобы использовать данный метод необходимо установить значение MPLEVL=2 в группе $CONTRL. Референтная конфигурация может быть получена методом RHF, ROHF или UHF (см. метод HF). Для данного метода доступно аналитическое вычисление градиента. Общая настройка метода производится в группе $MP2:

В GAMESS (US) метод MP2, по-умолчанию, не применяется по отношению к остовным орбиталям. Число остовных \(\alpha\)-орбиталей задаётся параметром NACORE, а число остовных \(\beta\)-орбиталей – параметром NBCORE.

Примечание

Некоторые свойства системы могут быть получены на основе градиента её энергии, однако его вычисление затратнее, чем расчёт энергии, поэтому в GAMESS (US) по-умолчанию они рассчитываются только в градиентных расчётах, так как это практически не влияет на общие затраты расчёта.

Чтобы включить расчёт свойств в рассчёт методом MP2, необходимо установить значение MP2PRP=.TRUE. в группе $MP2.

Метод MP2, в котором референтная конфигурация получена методом RHF, часто обозначается как RMP2. Аналогично существуют обозначения ROMP2 и UMP2.

В отличие от методов RMP2 и UMP2, для которых имеется чёткая методология, метод ROMP2 может быть реализован по-разному и, более того, различные его реализации приводят к немного разным результатам. В GAMESS (US) есть возможность задать его в трёх вариантах:

• RMP

С технической точки зрения данный вариант максимально близок к UMP2: выражения для \(\alpha\)- и \(\beta\)-орбиталей ищутся раздельно с помощью систем уравнений Попла-Несбета, а не системы Рутаана. Однако для занятых \(\alpha\)-орбиталей и соответствующих им \(\beta\)-орбиталей приравниваются пространственные части (\(\Psi^\alpha_i=\Psi^\beta_i\)), что и приводит к результату ROMP2.

Задаётся он следующим образом:

 $CONTRL
SCFTYP=ROHF MPLEVL=2
MULT=N
 $END
 $MP2
OSPT=RMP
 $END

Имеется также альтернативный способ задачи данного метода:

 $CONTRL
SCFTYP=UHF MPLEVL=2
MULT=N
 $END
 $SCF CUHF=.TRUE. $END

• ZAPT

Данный вариант, наоборот, ближе к RMP2: выражения для \(\alpha\)- и \(\beta\)-орбиталей ищутся совместно с помощью системы Рутаана при условии, что пространственные части соответствующих спин-орбиталей равны.

Задаётся он следующим образом:

 $CONTRL
SCFTYP=ROHF MPLEVL=2
MULT=N
 $END
 $MP2 OSPT=ZAPT $END

• OPT1

В этом варианте для соответствующих друг другу остовных \(\alpha\)- и \(\beta\)-орбиталей пространственные части приравниваются, но остовные и активные орбитали описываются раздельно.

Задать данный вариант ROMP2 можно следующим образом:

 $CONTRL
SCFTYP=MCSCF MPLEVL=2
MULT=N
 $END
 $DET
NCORE=M
NACT=N-1       !Все активные электроны – это неспаренные электроны
NELS=N-1       !Число активных орбиталей равно числу активных электронов
SZ=(N-1)/2     !Мультиплетность максимальна
PURES=.TRUE.   !И фиксирована
 $END          ! => Только одна конфигурация

Метод OPT1 лучше всего описывает систему (особенно, с точки зрения спина, ведь он в точности равен заданному), однако в GAMESS (US) его возможности крайне ограничены (например, невозможны градиентные расчёты).

Метод RMP – менее точный и довольно затратный (относительно ZAPT). Кроме того, в GAMESS (US), для недоступны градиентные и параллельные расчёты.

Метод ZAPT в GAMESS (US) поддерживает как градиентные так и параллельные расчёты, при этом по точности он не хуже RMP, но уступает OPT1.

Поэтому по умолчанию в GAMESS (US) устанолен метод ZAPT, то есть OSPT=ZAPT в группе $MP2.

В GAMESS (US) имеется шесть реализаций кода метода MP2, которые отличаются между собой по требуемым вычислительным ресурсам (дисковой и оперативной памяти):

• CODE=SERIAL

• CODE=DDI

• CODE=IMS

• CODE=RIMP2

• CODE=OMPRIMP2

• CODE=GPURIMP2

Примечание

Вышеперечисленные значения соответствуют параметру CODE в группе $MP2.

Важно

Фактически, первые три значения соответствуют методу $MP2, в то время как последние три соответствуют методу RI-MP2, который не только отличается требуемым объёмом вычислительных ресурсов, но и энергией, получаемой в результате расчёта.

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

	SCFTYP=RHF	SCFTYP=UHF	SCFTYP=ROHF OSPT=ZAPT	SCFTYP=ROHF OSPT=RMP
CODE=SERIAL	G	G	E	E
CODE=DDI	G	G	G	–
CODE=IMS	G	–	–	–
CODE=RIMP2	E	E	–	–
CODE=OMPRIMP2	E	–	–	–
CODE=GPURIMP2	E	–	–	–

Примечание

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

Теперь рассмотрим их несколько подробнее:

Первые три, как упоминалось выше, соответствуют методу MP2 в чистом виде. Отличаются они друг от дрга, в первую очередь объёмом требуемой оперативной и дисковой памяти.

Обозначим как \(N\) число базисных функций, а как \(M\) – число орбиталей, к которым применяется метод MP2 (ввиду наличия остовных орбиталей, к которым MP2 не применяется, \(M\) меньше, чем \(N\)).

Итак, при CODE=SERIAL требуемый объём оперативной памяти пропорционален \(N^3\), при этом расчёт проводится только последовательно и требует достаточно большой объём дисковой памяти. CODE=SERIAL является значением по умолчанию для последовательных расчётов при значениях SCFTYP=UHF или SCFTYP=ROHF.

При CODE=IMS требуемый объём оперативной памяти пропорционален \(N{\cdot}M^2\), при этом также требуется большой объём дисковой памяти, так как все данные в процессе расчёта сохраняются на локальные диски (однако для последовательных расчётах требования к диску ниже, чем при CODE=SERIAL). В целом, при данном значении последовательные расчёты идут быстрее, чем при CODE=SERIAL, помимо этого доступны и параллельные расчёты. CODE=IMS является значением по умолчанию при значении SCFTYP=RHF.

Стоит иметь ввиду, что запись и считывание с диска могут играть роль «бутылочного горлышка» в расчёте при большом объёме данных, хранимых на диске, то есть являтся причиной его замедления, ввиду простаивания вычислительных мощностей

При CODE=DDI требуемый объём оперативной памяти пропорционален \(N^4\), однако она распределяется между всеми узлами, кроме того, требования к дисковой памяти значительно ниже, а запись и считывание с диска практически отсутствуют. CODE=DDI является значением по умолчанию для параллельных расчётов при значениях SCFTYP=UHF или SCFTYP=ROHF.

Так как при CODE=DDI используется распределённая оперативная память, то требуется указать выделяемый её объём с помощью параметра MEMDDI в группе $SYSTEM (значения по умолчанию нет). Объём указывается в мегасловах (\(1\phantom{n}мегаслово=8{\cdot}10^6\phantom{n}байт\)).

Итак, последние три значения параметра CODE: CODE=RIMP2, CODE=OMPRIMP2 и CODE=GPURIMP2 соответствуют методу RI-MP2, при этом последние два являются, соответственно, реализациями метода RI-MP2 с использованием OpenMP и с вовлечением в расчёт графических процессоров.

Итак, в чём же заключается метод RI-MP2?

Для этого рассмотрим выражение энергии в методе MP2:

\[E^{\mathbf{MP2}}_0=E^{(0)}_0+E^{(1)}_i+E^{(2)}_0= \langle{0}|{\hat{H}_0}|{0}\rangle+\langle{0}|{\hat{V}}|{0}\rangle+ \sum_{j}'\frac{\langle{0}|{\hat{V}}|{D_i}\rangle\langle{D_i}|{\hat{V}}|{0}\rangle}{E^{(0)}_0-E^{(0)}_i}\]

Как видно из выражения выше, большая часть интегралов имеет вид:

\[\langle{0}|{\hat{V}}|{D_i}\rangle=\int\Psi^*_0\hat{V}\Psi^{rs}_{ab}dx\]

Именно вычисление данных интегралов и требует большую часть ресурсов.

Так как оператор возмущения \(\hat{V}\) может быть сведён к сумме двухэлектронных операторов, о чём упоминалось выше, то для него применимо одно из правил Слэтера-Кондона, согласно которому:

\[\begin{split}\int\Psi^*_0\hat{V}\Psi^{rs}_{ab}dx&= \\ =\int\int\psi^*_a(x_1)\psi^*_b(x_2)\hat{V}\psi_r(x_1)\psi_s(x_2)dx_1dx_2&- \int\int\psi^*_a(x_1)\psi^*_b(x_2)\hat{V}\psi_r(x_2)\psi_s(x_1)dx_1dx_2\end{split}\]

Примечание

Первый и второй интегралы отличаются друг от друга только перестановкой координат электронов в выражениях спин-орбиталей \(\psi_r\) и \(\psi_s\).

Таким образом, расчёт сводится к вычислению интегралов, зависящих от четырёх спин-орбиталей. Такие интегралы именуются четырёхцентровыми:

\[I_4=\int\int\psi^*_a(x_1)\psi^*_b(x_2)\hat{V}\psi_r(x_1)\psi_s(x_2)dx_1dx_2=\int\int\psi^*_a(x_1)\psi_r(x_1)\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx_1dx_2\]

Однако, это очень трудоёмко. Выходом из данной ситуации является применение техники, именуемой в англоязычной литературе Resolution of the Identity или, сокращённо, RI.

Рассмотрим принцип работы данной техники. Пусть существует набор функций \(\{\Theta_i\}\), такой что выполняется равенство:

\[\langle{\Theta}|\mathbf{\hat{V}}|\Theta\rangle=\sum_{ij}\langle{\Theta_i}|\hat{V}|\Theta_j\rangle=1\]

Запишем в развёрнутой форме данное выражение:

\[\langle{\Theta}|\mathbf{\hat{V}}|\Theta\rangle=\sum_{ij}\int\int\Theta^*_i(x_1)\hat{V}\Theta_j(x_2)dx_1dx_2\]

Представим его в матричной форме:

\[\begin{split}\sum_{ij}\int\int\Theta^*_i(x_1)\hat{V}\Theta_j(x_2)dx_1dx_2= \begin{pmatrix} \int\int\Theta^*_1(x_1) & \int\int\Theta^*_2(x_1) & \ldots \\ \end{pmatrix} \mathbf{\hat{V}} \begin{pmatrix} \Theta_1(x_2)dx_1dx_2 \\ \Theta_2(x_2)dx_1dx_2 \\ \vdots \\ \end{pmatrix}\end{split}\]

Таким образом, в данном выражении:

\[\begin{split}\langle{\Theta}|= \begin{pmatrix} \int\int\Theta^*_1(x_1) & \int\int\Theta^*_2(x_1) & \ldots \\ \end{pmatrix}\end{split}\]

\[\begin{split}|{\Theta}\rangle= \begin{pmatrix} \Theta_1(x_2)dx_1dx_2 \\ \Theta_2(x_2)dx_1dx_2 \\ \vdots \\ \end{pmatrix}\end{split}\]

\[\begin{split}\mathbf{\hat{V}}= \begin{pmatrix} \hat{V} & \hat{V} & \ldots \\ \hat{V} & \hat{V} & \ldots \\ \vdots & \vdots & \ddots \\ \end{pmatrix} =\hat{V} \begin{pmatrix} 1 & 1 & \ldots \\ 1 & 1 & \ldots \\ \vdots & \vdots & \ddots \\ \end{pmatrix} =\hat{V}\mathbf{I}\end{split}\]

Учитывая, что \(\hat{V}\mathbf{I}=\mathbf{I}\hat{V}\), получим:

\[\begin{split}\sum_{ij}\int\int\Theta^*_i(x_1)\hat{V}\Theta_j(x_2)dx_1dx_2= \begin{pmatrix} \int\int\Theta^*_1(x_1) & \int\int\Theta^*_2(x_1) & \ldots \\ \end{pmatrix} \mathbf{I}\hat{V} \begin{pmatrix} \Theta_1(x_2)dx_1dx_2 \\ \Theta_2(x_2)dx_1dx_2 \\ \vdots \\ \end{pmatrix}\end{split}\]

Теперь внесём оператор \(\hat{V}\) в столбец:

\[\begin{split}\sum_{ij}\int\int\Theta^*_i(x_1)\hat{V}\Theta_j(x_2)dx_1dx_2= \begin{pmatrix} \int\int\Theta^*_1(x_1) & \int\int\Theta^*_2(x_1) & \ldots \\ \end{pmatrix} \mathbf{I} \begin{pmatrix} \hat{V}\Theta_1(x_2)dx_1dx_2 \\ \hat{V}\Theta_2(x_2)dx_1dx_2 \\ \vdots \\ \end{pmatrix}\end{split}\]

Транспонируем полученное выражение, учитывая, что \(\mathbf{I}^T=\mathbf{I}\):

\[\begin{split}\langle{\Theta}|\mathbf{\hat{V}}|\Theta\rangle^T= \hat{V}|\Theta\rangle^T\mathbf{I}^T\langle{\Theta}|^T= \begin{pmatrix} \hat{V}\Theta_1(x_2)dx_1dx_2 & \hat{V}\Theta_2(x_2)dx_1dx_2 & \ldots \\ \end{pmatrix} \mathbf{I} \begin{pmatrix} \int\int\Theta^*_1(x_1) \\ \int\int\Theta^*_2(x_1) \\ \vdots \\ \end{pmatrix}\end{split}\]

Так как \(\langle{\Theta}|\mathbf{\hat{V}}|\Theta\rangle=1\), то и \(\langle{\Theta}|\mathbf{\hat{V}}|\Theta\rangle^T=1\), а значит, данное выражение можно «вклинить» в интеграл \(I_4\), не меняя его значение:

\[\begin{split}&\int\int\psi^*_a(x_1)\psi_r(x_1)\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx_1dx_2= \\ =&\int\int\psi^*_a(x_1)\psi_r(x_1)1\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx_1dx_2= \\ =&\int\int\psi^*_a(x_1)\psi_r(x_1)\langle{\Theta}|\mathbf{\hat{V}}|\Theta\rangle^T\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx_1dx_2= \\ =&\int\int\psi^*_a(x_1)\psi_r(x_1)\hat{V}|{\Theta}\rangle^T\mathbf{I}\langle{\Theta}|^T\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx_1dx_2 \\\end{split}\]

Для краткости заменим \(dx_1dx_2\) на \(dx\), а двойной интеграл заменим одиночным:

\[I_4=\int\psi^*_a(x_1)\psi_r(x_1)\hat{V}|{\Theta}\rangle^T\mathbf{I}\langle{\Theta}|^T\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx\]

А теперь развернём выражение \(\hat{V}|{\Theta}\rangle^T\langle{\Theta}|^T\) в матричную форму:

\[\begin{split}\int\psi^*_a(x_1)\psi_r(x_1) \hat{V} \begin{pmatrix} \Theta_1(x_2)dx & \Theta_2(x_2)dx & \ldots \\ \end{pmatrix} \mathbf{I} \begin{pmatrix} \int\Theta^*_1(x_1) \\ \int\Theta^*_2(x_1) \\ \vdots \\ \end{pmatrix} \hat{V}\psi^*_b(x_2)\psi_s(x_2)dx\end{split}\]

А теперь внесём в строку всё, что стоит левее её, а в столбец – всё что правее его:

\[\begin{split}\begin{pmatrix} \int\psi^*_a(x_1)\psi_r(x_1)\hat{V}\Theta_1(x_2)dx & \int\psi^*_a(x_1)\psi_r(x_1)\hat{V}\Theta_2(x_2)dx & \ldots \\ \end{pmatrix} \mathbf{I} \begin{pmatrix} \int\Theta^*_1(x_1)\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx \\ \int\Theta^*_2(x_1)\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx \\ \vdots \\ \end{pmatrix}\end{split}\]

Для ещё большей краткости примем обозначения:

\[(ar|\Theta_i):=\int\psi^*_a(x_1)\psi_r(x_1)\hat{V}\Theta_i(x_2)dx\]

\[(\Theta_i|bs):=\int\Theta^*_i(x_1)\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx\]

В итоге получаем:

\[\begin{split}I_4= \begin{pmatrix} (ar|\Theta_1) & (ar|\Theta_2) & \ldots \\ \end{pmatrix} \mathbf{I} \begin{pmatrix} (\Theta_1|bs) \\ (\Theta_2|bs) \\ \vdots \\ \end{pmatrix}\end{split}\]

Или, после перемножения:

\[I_4=\sum_{ij}(ar|\Theta_i)(\Theta_j|bs)\]

Что в развёрнутой форме равно:

\[I_4=\sum_{ij}\int\int\psi^*_a(x_1)\psi_r(x_1)\hat{V}\Theta_i(x_2)dx_1dx_2\int\int\Theta^*_j(x_1)\hat{V}\psi^*_b(x_2)\psi_s(x_2)dx_1dx_2\]

Таким образом, применив технику RI, можно «разбить» интеграл, зависящий от четырёх спин-орбиталей (четырёхцентровый \(I_4\)) на сумму произведений интегралов, зависящих только от двух спин-орбиталей и одной функции \(\Theta_i\) (трёхцентровых \(I_3\)). Подобный приём позволяет существенно снизить вычислительную стоимость расчёта.

Единственной сложностью в реализации данной техники является наличие таких функций \(\Theta_i\), для которых выполняется условие:

\[\sum_{ij}\langle{\Theta_i}|\hat{V}|\Theta_j\rangle=1\]

Поэтому на практике вместо набора \(\{\Theta_i\}\) используется набор взаимно ортогональных функций \(\{\Delta_i\}\), а несоблюдение указанного выше условия компенсируется делением каждого слагаемого на \(\langle{\Delta_i}|\hat{V}|\Delta_j\rangle\). Так как \(\langle{\Delta_i}|\hat{V}|\Delta_j\rangle\) зависит только от двух функций (\(I_2\)), то его наличие в итоговом практически не влияет на объём требуемых вычислительных ресурсов.

Таким образом:

\[I_4\approx\sum_{ij}\frac{(ab|\Delta_i)(\Delta_j|rs)}{\langle{\Delta_i}|\hat{V}|\Delta_j\rangle}\]

Стоит отметить, что данное выражение является лишь приближением к истинному результату, поэтому расчёты с использованием техники RI чувствительны к набору функций \(\{\Delta_i\}\).

Данный набор по аналогии со молекулярными орбиталями может быть разложен в виде линейной комбинации функций, привязанных к каждому атому, – аналогов атомных орбиталей. Эти функции, как и атомные орбитали, могут быть представлены в виде разложения по базисному набору \(\{\delta_i\}\), привязанному к конкретному сорту атомов.

Набор функций \(\{\delta_i\}\) в англоязычной литературе именуется как Auxilary Basis Set, то есть вспомогательный базисный набор. Как правило, вспомогательные базисные наборы разрабатываются под какие-то конкретные обыкновенные базисные наборы.

В GAMESS (US) вспомогательный базисный набор задаётся в группе $AUXBAS, при этом имеется два способа:

Первый способ – это с помощью параметра CABNAM: он позволяет выбрать небольшое число вспомогательных базисов. Доступны следующие вспомогательные базисные наборы:

• CABNAM=SVP (def2-SVP-RIFIT) – оптимизирован под def2-SVP.

• CABNAM=TZVP (def2-TZVP-RIFIT) – оптимизирован под def2-TZVP.

• CABNAM=TZVPP (def2-TZVPP-RIFIT) – оптимизирован под def2-TZVPP.

• CABNAM=CCD (cc-pVDZ-RIFIT) – оптимизирован под cc-pVDZ.

• CABNAM=ACCD (aug-cc-pVDZ-RIFIT) – оптимизирован под aug-cc-pVDZ.

• CABNAM=CCT (cc-pVTZ-RIFIT) – оптимизирован под cc-pVTZ.

• CABNAM=ACCT (aug-cc-pVTZ-RIFIT) – оптимизирован под aug-cc-pVTZ.

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

Второй способ – это установить значение EXTCAB=.TRUE. и с помощью параметра CABNAM указать название базисного набора, записанного во внешнем файле. Формат записи вспомогательных базисных наборов во внешнем файле не отличается от такового для обычных базисов.

По-умолчанию, в GAMESS (US) расчёт методом RI-MP2 протекает с использованием распределённой оперативной памяти. А значит, для него необходима установка выделяемого её объёма (с помощью параметра MEMDDI в группе $SYSTEM). При желании можно с помощью значений GOSMP=.FALSE. и USEDM=.FALSE. в группе $RIMP2 перейти к локализованной оперативной памяти.

Важно

Параллельный расчёт методом RI-MP2 использует также особый тип памяти, который не регулируется параметрами. Чтобы удостовериться в достаточности объёма выделенной памяти, рекомендуется предварительно запустить задачу со значением EXETYP=CHECK в группе $CONTRL.