Метод самосогласованного поля

Метод самосогласованного поля (SCF) – это основополагающий метод квантовой химии, который позволяет на основе данных о составе системы (\(q^{nuc}_i\), \(Z^{nuc}_i\), \(Z_\sum\), \(s^{elec}_i\)), а также базисного набора \(\{\phi_p\}\) составить её приближённое описание, из которого можно извлечь её свойтства системы. Частными случаями метода SCF являются:

• Метод Хартри-Фока (HF) – метод SCF на основе волновой функции системы \(\Psi\).

• Метод Кона-Шэма (KS) – метод SCF на основе электронной плотности системы \(\rho\).

На этих двух методах основаны практически все остальные квантово-химические методы.

Метод Хартри-Фока

Метод Хартри-Фока позволяет построить аппроксимацию волновой функции системы \(\Psi\) с помощью базисного набора \(\{\phi_p\}\).

Перед тем как перейти к краткому рассмотрению метода HF рассмотрим стационарное уравнение Шрёдингера:

\[\hat{\mathcal{H}}\Psi=E\Psi\]

В нём присутствует полная энергия системы \(E\). Но из чего она состоит?

Полная энергия системы складывается из потенциальных энергий взаимодействия всех частиц в её составе \(U\), то есть, электронов и ядер, а также их кинетических энергий \(T\). Поэтому её можно представить в виде суммы пяти слагаемых:

\[E=T^{n}+T^{e}+U^{nn}+U^{ne}+U^{ee}\]

где:

• \(T^{n}\) – кинетическая энергия ядер.

• \(T^{e}\) – кинетическая энергия электронов.

• \(U^{nn}\) – потенциальная энергия взаимодействия ядер друг с другом.

• \(U^{ne}\) – потенциальная энергия взаимодействия ядер и электронов.

• \(U^{ee}\) – потенциальная энергия взаимодействия электронов друг с другом.

Абсолютно аналогично можно разделить и гамильтониан \(\hat{\mathcal{H}}\). Собственными значениями полученных операторов и будут вышеприведённые состовляющие полной энергии:

\[\hat{\mathcal{H}}=\hat{T}^{n}+\hat{T}^{e}+\hat{U}^{nn}+\hat{U}^{ne}+\hat{U}^{ee}\]

Взаимная зависимость координат электронов и ядер делает задачу нахождения волновой функции \(\Psi\) практически неразрешимой, так как это исключает возможность разделения переменных (координат электронов и ядер).

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

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

\[ \begin{align}\begin{aligned}\hat{\mathcal{H}}&=\underbrace{\hat{T}^{n}+\hat{U}^{nn}}+\underbrace{\hat{T}^{e}+\hat{U}^{ne}}+\underbrace{\phantom{|}\hat{U}^{ee}\phantom{|}}\\\hat{\mathcal{H}}&=\phantom{)))}\hat{H}{(0e)}\phantom{)}+\phantom{M}\hat{H}{(1e)}\phantom{|}+\hat{H}{(ne)}\end{aligned}\end{align} \]

Таким образом, гамильтониан состоит из трёх частей:

• \(\hat{H}{(0e)}\) – независящая от расположения электронов.

• \(\hat{H}{(1e)}\) – определяющаяся расположением каждого электрона в отдельности.

• \(\hat{H}{(ne)}\) – определяющаяся взаимным расположением электронов друг относительно друга.

Для решения уравнения Шрёдингера необходимо разделить переменные (координаты электронов):

Для первого слагаемого эта задача тривиальна, так как оно не зависит от них.

Второе слагаемое можно представить в виде суммы одноэлектронных операторов \(\hat{h}(q_i)\)

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

Основная идея метода Хартри-Фока – это замена слагаемого \(\hat{H}{(ne)}\) на сумму потенциалов электронов \(\hat{v}(q_i)\) в усреднённом поле, создаваемом другими электронами. Данные потенциалы именуются потенциалами Хартри-Фока.

Примечание

Потенциал Хартри-Фока может быть разбит на две группы слагаемых:

\[\hat{v}(x_i)=\sum_{j}\hat{J}_j(q_i)-\sum_{j}\hat{K}_j(q_i)\]

где:

• \(\hat{J}_j(x_i)\) – кулоновский оператор, собственным значением которого является энергия отталкивания электрона \(i\), расположенного на спин-орбитали \(\Psi_j\), от остальных электронов.

• \(\hat{K}_j(x_i)\) – обменный оператор, собственным значением которого является понижение энергии вследствие представления волновой функции в виде определителя Слэтера. То есть, он имеет чисто математическую природу.

Таким образом, гамильтониан в методе Хартри-Фока (\(\hat{\mathcal{H}}^{\mathbf{HF}}\)) имеет вид:

\[\hat{\mathcal{H}}^{\mathbf{HF}}=\hat{H}{(0e)}+\sum_{i}\hat{h}(q_i)+\sum_{i}\hat{v}(q_i)\]

Или:

\[\hat{\mathcal{H}}^{\mathbf{HF}}=\hat{H}{(0e)}+\sum_{i}\hat{f}(q_i)\]

где \(\hat{f}(q_i)=\hat{h}(q_i)+\hat{v}(q_i)\) – одноэлектронный оператор Фока.

Примечание

Таким образом, собственными значениями одноэлектронного оператора Фока \(\hat{f}(q_i)\) являются орбитальные энергии \(\varepsilon_j\). Орбитальная энергия \(\varepsilon_j\) – это сумма кинетической энергии электрона, находящегося на спин-орбитали \(\psi_j\), энергии его взаимодействия с ядрами атомов и энергии его взаимодействия с усреднённым полем, создаваемым другими электронами.

Так как слагаемое \(\hat{H}{(0e)}\) не зависит от координат электронов, то оно может быть рассчитано отдельно от остальной части, поэтому обычно оперируют не полным гамильтонианом \(\hat{\mathcal{H}}^{\mathbf{HF}}\), а электронным гамильтонианом \(\hat{H}^{\mathbf{HF}}\):

\[\hat{H}^{\mathbf{HF}}=\sum_{i}\hat{f}(q_i)\]

Примечание

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

Итого: метод Хартри-Фока позволяет избавиться от взаимной зависимости координат электронов, что позволяет найти решение уравнения Шрёдингера в приближённом виде.

Более подробно, метод Хартри-Фока работает следующим образом:

1) Выражение для атомных орбиталей \(\chi_j\) ищется в виде линейной комбинации базисных функций \(\phi_p\):

\[\chi_j=\sum^P_{p=1}d_{jp}\phi_p\]

Коэффициенты \(d_{jp}\) находятся вариционным методом Ритца так, чтобы энергия атомной орбитали \(\chi_j\) была минимальной.

Примечание

В общем виде, вариационный метод Ритца заключается в поиске минимума энергии \(E\) как функции от коэффициентов \(c_i\) перед базисными функциями \(\varphi_i\) в разложении волновой функции \(\Psi\):

\[\Psi=\sum^n_{i=1}c_i\varphi_i\]

Волновая функция \(\Psi\) является собственной функцией гамильтониана \(\hat{H}\), соответствующей собственному значению \(E\). Таким образом, справедливо уравнение Шрёдингера следующего вида:

\[\hat{H}\Psi=E\Psi\]

При подстановке в него разложения волновой функции \(\Psi\) получаем выражение:

\[\hat{H}\sum^n_{i=1}c_i\varphi_i=E\sum^n_{i=1}c_i\varphi_i\]

Умножим данное выражение на разложение волновой функции \(\Psi^*\), комплексно сопряжённой к \(\Psi\), и проинтегрируем по совокупности переменных \(x\), от которых зависят данные волновые функции:

\[\int\sum^n_{j=1}c^*_j\varphi^*_j\hat{H}\sum^n_{i=1}c_i\varphi_idx=\int\sum^n_{j=1}c^*_j\varphi^*_j\sum^n_{i=1}c_i\varphi_idx\]

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

\[\int\sum^n_{j=1}\sum^n_{i=1}c^*_jc_i\varphi^*_j\hat{H}\varphi_idx=\int{E}\sum^n_{i=1}\sum^n_{j=1}c^*_jc_i\varphi^*_j\varphi_idx\]

\[\sum^n_{j=1}\sum^n_{i=1}c^*_jc_i\int\varphi^*_j\hat{H}\varphi_idx=E\sum^n_{j=1}\sum^n_{i=1}c^*_jc_i\int\varphi^*_j\varphi_idx\]

Продифференцируем полученное выше выражение по \(c_k\):

\[\frac{\partial}{\partial{c_k}}\sum^n_{j=1}\sum^n_{i=1}c^*_jc_i\int\varphi^*_j\hat{H}\varphi_idx= \frac{\partial}{\partial{c_k}}\left(E\sum^n_{j=1}\sum^n_{i=1}c^*_jc_i\int\varphi^*_j\varphi_idx\right)\]

\[\sum^n_{j=1}c^*_j\int\varphi^*_j\hat{H}\varphi_idx= \frac{\partial{E}}{\partial{c_k}}\sum^n_{j=1}\sum^n_{i=1}c^*_jc_i\int\varphi^*_j\varphi_idx+ E\sum^n_{j=1}c^*_j\int\varphi^*_j\varphi_idx\]

Необходимым условием минимума энергии является равенство нулю всех первых производных энергии \(E\) по коэффициентам \(c_k\):

\[\frac{\partial{E}}{\partial{c_k}}=0,{\forall}k\in\{1,2,3,\ldots,n\}\]

Тогда:

\[\sum^n_{j=1}c^*_j\int\varphi^*_j\hat{H}\varphi_idx=E\sum^n_{j=1}c^*_j\int\varphi^*_j\varphi_idx\]

\[\sum^n_{j=1}\left(c_j\int\varphi^*_i\hat{H}\varphi_jdx\right)^*=E\sum^n_{j=1}\left(c_j\int\varphi^*_i\varphi_jdx\right)^*\]

\[\sum^n_{j=1}c_j\int\varphi^*_i\hat{H}\varphi_jdx=E\sum^n_{j=1}c_j\int\varphi^*_i\varphi_jdx\]

Перенеся всё в левую часть, получаем систему уравнений:

\[\sum^n_{j=1}c_j\int\varphi^*_i\hat{H}\varphi_jdx-E\sum^n_{j=1}c_j\int\varphi^*_i\varphi_jdx=0,{\forall}i\in\{1,2,3,\ldots,n\}\]

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

\[\begin{split}\begin{pmatrix} H_{11}-ES_{11} & H_{12}-ES_{12} & \ldots & H_{1n}-ES_{1n} \\ H_{21}-ES_{21} & H_{22}-ES_{22} & \ldots & H_{2n}-ES_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ H_{n1}-ES_{n1} & H_{n2}-ES_{n2} & \ldots & H_{nn}-ES_{nn} \\ \end{pmatrix} \times\ \begin{pmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \\ \end{pmatrix} =\mathbf{0}\end{split}\]

где:

• \(H_{ij}=\int\varphi^*_i\hat{H}\varphi_jdx\) – матричный элемент гамильтониана.

• \(S_{ij}=\int\varphi^*_i\varphi_jdx\) – интеграл перекрывания.

• \(\mathbf{0}\) – нулевой столбец, высотой \(n\).

Данная система уравнений является однородной, а значит, имеет тривиальное решение \(c_i=0\). Существование нетривиальных решений для данной системы уравнение возможно, только если выполняется следующее равенство:

\[\begin{vmatrix}H_{ij}-ES_{ij}\end{vmatrix}=0\]

которое в развёрнутом виде выглядит следующим образом:

\[\begin{split}\begin{vmatrix} H_{11}-ES_{11} & H_{12}-ES_{12} & \ldots & H_{1n}-ES_{1n} \\ H_{21}-ES_{21} & H_{22}-ES_{22} & \ldots & H_{2n}-ES_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ H_{n1}-ES_{n1} & H_{n2}-ES_{n2} & \ldots & H_{nn}-ES_{nn} \\ \end{vmatrix} =0\end{split}\]

Данное равенство также известно как вековое уравнение. Вековое уравнение имеет \(n\) решений, наименьшее из которых – это энергия основного состояния \(E^0\), а остальные \(E^{k\ne0}\) – энергии возбуждённых.

Подстановка \(E^k\) в приведённую выше систему уравнений даёт набор коэффициентов \(\{c^k_i\}\) в разложении волновой функции \(\Psi^k\).

2) Молекулярные орбитали \(\Psi_i\) представляются в виде линейной комбинации атомных орбиталей \(\chi_j\):

\[\Psi_i=\sum^J_{j=1}c_{ij}\chi_j\]

3) Каждая молекулярная спин-орбиталь \(\psi_i\) представляется в виде произведения пространственной части (то есть, орбитали \(\Psi_i\)) и спиновой части \(\varsigma_i\):

\[ \begin{align}\begin{aligned}\psi_{2i-1}=\Psi_i\alpha\\\psi_{2i}=\Psi_i\beta\end{aligned}\end{align} \]

где:

• \(\alpha\) – спиновая часть спин-орбитали, соответствующей спину, равному \(+\frac{1}{2}\hbar\) (\(\alpha\)-орбитали).

• \(\beta\) – спиновая часть спин-орбитали, соответствующей спину, равному \(-\frac{1}{2}\hbar\) (\(\beta\)-орбитали).

• \(\Psi_i\) – пространственная часть спин-орбитали.

4) Волновая функция системы \(\Psi\) представляется через молекулярные спин-орбитали \(\psi_k\) в виде определителя Слэтера (SD):

\[\begin{split}\Psi=\frac{1}{\sqrt{n!}} \begin{vmatrix} \psi_1(q_1,\sigma_1) & \psi_1(q_2,\sigma_2) & \ldots & \psi_1(q_n,\sigma_n) \\ \psi_2(q_1,\sigma_1) & \psi_2(q_2,\sigma_2) & \ldots & \psi_2(q_n,\sigma_n) \\ \vdots & \vdots & \ddots & \vdots \\ \psi_n(q_1,\sigma_1) & \psi_n(q_2,\sigma_2) & \ldots & \psi_n(q_n,\sigma_n) \\ \end{vmatrix}\end{split}\]

Примечание

\(q_j\) – это пространственные координаты электрона под номером \(j\), а \(\sigma_j\) – это соответствующая ему спиновая переменная. У неё нет, как такового, физического смысла, но в данном случае это и не важно: главное, что пространственная часть спин-орбитали \(\Psi_i\) – это функция координат \(q_j\), а спиновая часть \(\varsigma_i\) – функция спиновой переменной \(\sigma_j\).

5) Для получения выражения волновой функции системы \(\Psi\) остаётся только определить коэффициенты \(c_{ij}\) в разложении молекулярных орбиталей.

Рассмотрим одноэлектронный оператор Фока \(\hat{f}\):

\[\hat{f}(q_i)=\hat{h}(q_i)+\hat{v}(q_i)\]

Так как собственными значениями оператора Фока являются орбитальные энергии \(\varepsilon_j\) электронов на орбиталях \(\psi_j\), которые являются его собственными функциями, то:

\[\hat{f}(q_i)\psi_j(q_i,\sigma_i)=\varepsilon_j\psi_j(q_i,\sigma_i)\]

Следовательно:

\[\varepsilon_i=\int\int\psi^*_i(q,\sigma)\hat{f}(q)\psi_i(q,\sigma)dqd\sigma\]

Для нахождения коэффициентов \(c_{ij}\) составляется система уравнений Рутаана:

\[\mathbf{F}\times\mathbf{C}=\mathbf{E}\times\mathbf{S}\times\mathbf{C}\]

где:

• \(\mathbf{F}\) – матрица Фока, с элементами \(F_{ij}=\int\int\psi^*_i(q,\sigma)\hat{f}(q)\psi_j(q,\sigma)dqd\sigma\).

• \(\mathbf{S}\) – матрица перекрывания, с элементами \(S_{ij}=\int\int\psi^*_i(q,\sigma)\psi_j(q,\sigma)dqd\sigma\).

• \(\mathbf{C}\) – матрица коэффициентов, с элементами \(C_{ij}=c_{ij}\).

• \(\mathbf{E}\) – матрица орбитальных энергий, с элементами \(E_{ij}=\varepsilon_i\) при \(i=j\), \(E_{ij}=0\) при \(i\ne{j}\).

Примечание

Далее для сокращения записи (там, где это не принципиально) пространственная \(q\) и спиновая \(\sigma\) переменные будут объединены в обобщённую переменную \(x\). По этой же причине знак двойного интегрирования будет заменён одиночным, например:

\[\int\Psi^*_i\hat{H}\Psi_jdx:=\int\int\Psi^*_i\hat{H}\Psi_jdqd\sigma\]

В развёрнутом виде данная система уравнений выглядит следующим образом:

\[\begin{split} \begin{pmatrix} F_{11}-\varepsilon_1S_{11} & F_{12}-\varepsilon_1S_{12} & \ldots & F_{1n}-\varepsilon_1S_{1n} \\ F_{21}-\varepsilon_2S_{21} & F_{22}-\varepsilon_2S_{22} & \ldots & F_{2n}-\varepsilon_2S_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ F_{n1}-\varepsilon_nS_{n1} & F_{n2}-\varepsilon_nS_{n2} & \ldots & F_{nn}-\varepsilon_nS_{nn} \\ \end{pmatrix} \times \begin{pmatrix} c_{11} & c_{12} & \ldots & c_{1n} \\ c_{21} & c_{22} & \ldots & c_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ c_{n1} & c_{n2} & \ldots & c_{nn} \\ \end{pmatrix} =\mathbf{0}\end{split}\]

где \(\mathbf{0}\) – нулевая матрица размера \(n{\times}n\).

Итого: имеется система из \(n^2\) уравнений и \(n^2\) неизвестных \(c_{ij}\), из которой можно найти искомые коэффициенты.

Но для этого требуется определить для электрона \(i\) усреднённое поле, создаваемое другими электронами. Поле, создаваемое электронами – это не что иное, как электронная плотность системы \(\rho\):

\[\rho=\sum^n_{i=1}\Psi^*_i\Psi_i\]

Но, если разложить молекулярную орбиталь \(\Psi_i\) по базису из атомных орбиталей:

\[\rho=\sum^n_{i=1}\Psi^*_i\Psi_i=\sum^n_{i=1}\sum^n_{j=1}c^*_{ji}\chi^*_j\sum^n_{k=1}c_{ki}\chi_k= \sum^n_{j=1}\sum^n_{k=1}\sum^n_{i=1}c^*_{ji}c_{ki}\chi^*_j\chi_k\]

то окажется, что электронная плотность зависит от искомых коэффициентов \(c_{ij}\).

Это деляет невозможным их точное нахождение, однако позволяет итерационно найти приближённые значения коэффициентов \(c_{ij}\), так как электронную плотность связывает с потенциалом Хартри-Фока матрица плотности с элементами \(P_{jk}\):

\[P_{jk}=\sum^n_{i=1}c^*_{ji}c_{ki}\]

Итак, чтобы найти приближённое значений коэффициентов \(c_{ij}\), необходимо:

• Приняв довольно грубые упрощения, получить набор молекулярных орбиталей с коэффициентами \(c^{0}_{ij}\).

• С помощью коэффициентов \(c^{0}_{ij}\) выразить элементы матрицы плотности.

• Зная выражение элементов матрицы плотности, решить систему уравнений Рутаана и получить новый набор коэффициентов \(c^{1}_{ij}\).

• С помощью коэффициентов \(c^{1}_{ij}\) снова выразить элементы матрицы плотности.

• Зная новое выражение элементов матрицы плотности, решить систему уравнений Рутаана и получить новый набор коэффициентов \(c^{2}_{ij}\).

• Продолжать до тех пор, пока изменение энергии всей системы при переходе от одного цикла к другому не опустится ниже некоторого порогового значения (иными словами, до схождения расчёта).

• Полученные на последней стадии коэффициенты \(c^{k}_{ij}\) и есть искомые коэффициенты в аппроксимации к волновой функции системы \(\Psi\).

Примечание

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

Разновидности метода Хартри-Фока

Строго говоря, спин-орбитали имеют различные не только спиновые части \(\varsigma\), но и пространственные части \(\Psi^\varsigma_i\). Поэтому более корректным будет записать выражения для спин-орбиталей в следующим виде:

\[ \begin{align}\begin{aligned}\psi_{2i-1}=\Psi^\alpha_i\alpha\\\psi_{2i}=\Psi^\beta_i\beta\end{aligned}\end{align} \]

где:

• \(\Psi^\alpha_i\) – пространственная часть \(\alpha\)-орбитали (спин равен \(+\frac{1}{2}\hbar\)).

• \(\Psi^\beta_i\) – пространственная часть \(\beta\)-орбитали (спин равен \(-\frac{1}{2}\hbar\)).

Неравенство \(\Psi^\alpha_i\ne\Psi^\beta_i\) делает расчёт более сложным, по сравнению с описанным выше случаем. В частности, система уравнений Рутаана разбивается на две системы уравнений Попла-Несбета (отдельные для \(\alpha\)- и \(\beta\)-орбиталей).

Однако приближение \(\Psi^\alpha_i=\Psi^\beta_i\equiv\Psi_i\) во многих случаях не вносит большой погрешности. В связи с этим, существует три основных разновидности метода Хартри-Фока:

Ограниченный метод Хартри-Фока (RHF) – данный метод основан на приближении \(\Psi^\alpha_i=\Psi^\beta_i\equiv\Psi_i\), которое позволяет существенно сократить затратность расчёта. Данный метод может использоваться только для описания систем с мультиплетностью, равной единице. То есть, только если все электроны в системе спарены.

Таким образом, определитель Слэтера для системы, содержащей \(2n\) электронов, в методе RHF принимает вид:

\[\begin{split}\Psi=\frac{1}{\sqrt{(2n)!}} \begin{vmatrix} \Psi_1(q_1)\alpha(\sigma_1) & \Psi_1(q_2)\alpha(\sigma_2) & \ldots & \Psi_1(q_{2n})\alpha(\sigma_{2n}) \\ \Psi_1(q_1)\beta(\sigma_1) & \Psi_1(q_2)\beta(\sigma_2) & \ldots & \Psi_1(q_{2n})\beta(\sigma_{2n}) \\ \vdots & \vdots & \ddots & \vdots \\ \Psi_n(q_1)\beta(\sigma_1) & \Psi_n(q_2)\beta(\sigma_2) & \ldots & \Psi_n(q_{2n})\beta(\sigma_{2n}) \\ \end{vmatrix}\end{split}\]

Ограниченный метод Хартри-Фока для систем с открытой оболочкой (ROHF) – данный метод позволяет описывать системы с неспаренными электронами. Спаренные электроны описываются таким же образом, как и в RHF. Неспаренные электроны, всегда располагаются на \(\alpha\)-орбиталях (спин равен \(+\frac{1}{2}\hbar\)), так как они всегда ниже по энергии, чем соответствующие \(\beta\)-орбитали.

Тогда определитель Слэтера для системы, содержащей \(2n\) спаренных электронов и \(m\) неспаренных электронов в методе ROHF принимает вид:

\[\begin{split}\Psi=\frac{1}{\sqrt{(2n+m)!}} \begin{vmatrix} \mathbf{M}^{2n}_{2n} & \mathbf{F}^{2n}_1 & \mathbf{F}^{2n}_2 & \ldots & \mathbf{F}^{2n}_m \\ \mathbf{A}^1_{2n} & \alpha^1_1 & \alpha^1_2 & \ldots & \alpha^1_m \\ \mathbf{A}^2_{2n} & \alpha^2_1 & \alpha^2_2 & \ldots & \alpha^2_m \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{A}^m_{2n} & \alpha^m_1 & \alpha^m_2 & \ldots & \alpha^m_m \\ \end{vmatrix}\end{split}\]

где приняты следующие обозначения:

\[\begin{split}\mathbf{M}^{2n}_{2n}:= \begin{pmatrix} \Psi_1(q_1)\alpha(\sigma_1) & \Psi_1(q_2)\alpha(\sigma_2) & \ldots & \Psi_1(q_{2n})\alpha(\sigma_{2n}) \\ \Psi_1(q_1)\beta(\sigma_1) & \Psi_1(q_2)\beta(\sigma_2) & \ldots & \Psi_1(q_{2n})\beta(\sigma_{2n}) \\ \vdots & \vdots & \ddots & \vdots \\ \Psi_n(q_1)\beta(\sigma_1) & \Psi_n(q_2)\beta(\sigma_2) & \ldots & \Psi_n(q_{2n})\beta(\sigma_{2n}) \\ \end{pmatrix}\end{split}\]

\[\begin{split}\mathbf{A}^i_{2n}:= \begin{pmatrix} \Psi_{2n+i}(q_1)\alpha(\sigma_1) & \Psi_{2n+i}(q_2)\alpha(\sigma_2) & \ldots & \Psi_{2n+i}(q_{2n})\alpha(\sigma_{2n}) \\ \end{pmatrix}\end{split}\]

\[\begin{split}\mathbf{F}^{2n}_j:= \begin{pmatrix} \Psi_1(q_{2n+j})\alpha(\sigma_{2n+j}) \\ \Psi_2(q_{2n+j})\beta(\sigma_{2n+j}) \\ \vdots \\ \Psi_{2n}(q_{2n+j})\beta(\sigma_{2n+j}) \\ \end{pmatrix}\end{split}\]

\[\alpha^i_j:=\Psi_{2n+i}(q_{2n+j})\alpha(\sigma_{2n+j})\]

Неограниченный метод Хартри-Фока (UHF), в отличие от RHF и ROHF не содержит приближение \(\Psi^\alpha_i=\Psi^\beta_i\), поэтому он более корректно описывает систему, но при этом требует большее количество вычислительных ресурсов. Как и ROHF, UHF позволяет описывать системы с неспаренными электронами.

В общем виде, определитель Слэтера для системы, содержащей \(n\) электронов, распределённых по \(m\) орбиталям, в методе UHF принимает вид:

\[\begin{split}\Psi=\frac{1}{\sqrt{n!}} \begin{vmatrix} \Psi^\alpha_1(q_1)\alpha(\sigma_1) & \Psi^\alpha_1(q_2)\alpha(\sigma_2) & \ldots & \Psi^\alpha_1(q_{n})\alpha(\sigma_{n}) \\ \Psi^\beta_1(q_1)\beta(\sigma_1) & \Psi^\beta_1(q_2)\beta(\sigma_2) & \ldots & \Psi^\beta_1(q_{n})\beta(\sigma_{n}) \\ \vdots & \vdots & \ddots & \vdots \\ \Psi^\varsigma_m(q_1)\varsigma(\sigma_1) & \Psi^\varsigma_m(q_2)\varsigma(\sigma_2) & \ldots & \Psi^\varsigma_m(q_{n})\varsigma(\sigma_{n}) \\ \end{vmatrix}\end{split}\]

Недостатки метода Хартри-Фока

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

Почему же это невозможно? Существует несколько причин:

Во-первых, как было описано выше, волновая функция системы \(\Psi\) сводится к волновым функциям ядер \(\Psi^{nuc}_i(q^{nuc})\) и электронов \(\Psi^{elec}_i(q^{elec})\), однако координаты ядер \(q^{nuc}\) и электронов \(q^{elec}\) зависят друг от друга, поэтому разделение этих переменных невозможно.

На помощь приходит приближение Борна-Оппенгеймера: так как массы ядер значительно больше масс электронов, то скорости движения ядер крайне малы, а значит ими можно пренебречь. Тогда координаты ядер \(q^{nuc}\) становятся постоянными и проблема разделения переменных исчезает.

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

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

И в-третьих, недостатки самого метода Хартри-Фока:

Настоящая квантово-химическая система не находится только в одном состоянии. Реальное её состояние является смесью всех состояний с преобладанием наиболее низких по энергии.

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

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

Данное явление именуется статической корреляцией.

Также в методе Хартри-Фока взаимодействия между электронами с параллельными спинами и взаимодействия между электронами с антипараллельными спинами описываются одинаково, что, вообще говоря, некорректно, ведь электроны с параллельными спинами не могут оказаться на одной орбитали \(\Psi_i\), а значит, они, в среднем, находятся на большем расстоянии друг от друга.

Данное явление именуется динамической корреляцией.

Однако оба вида корреляции могут быть учтемы теми или иными методами, многие из которых являются дальнейшим развитием метода Хартри-Фока.

Многоконфигурационные методы

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

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

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

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

Генерация возбуждённых конфигураций может быть реализована по-разному. Наиболее простой способ – это представление возбуждённых конфигураций в виде определителей Слэтера (SD), в которых строки с занятыми спин-орбиталями \(\psi_k\), заменены на строки с вакантными (в основном состоянии) орбиталями \(\psi_l\).

Альтернативный подход – это использование линейной комбинации данных определителей, также известной как функция конфигурационного состояния (CSF).

Таким образом, по количеству конфигураций данные методы могут быть разделены на:

• Одноконфигурационные [single-configurational] (SC).

• Многоконфигурационные [multiconfigurational] (MC).

А по количеству референтных конфигураций на:

• Однореферентные [single-reference] (SR).

• Многореферентные [multireference] (MR).

Разница между этими методами продемонстрирована на рисунке ниже:

Примечание

Все многореферентные методы одновременно являются многоконфигурационными, а все одноконфигурационные – однореферентными.

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

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

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

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

Таким образом, многоконфигурационные однореферентные методы учитывают динамическую корреляцию, но не учитывают статическую корреляцию.

Однако, если используемые конфигурации подвергнуть оптимизации по энергии, например, методом SCF, то ситуация изменится и вклады конфигураций будут соответствовать действительности. Именно на этом принципе основан, так называемый, многоконфигурационный метод самосогласованного поля (MCSCF).

Но, так как оптимизация каждой конфигурации – это, фактически, проведение всей процедуры метода SCF, то затратность метода MCSCF намного выше, чем метода HF, поэтому на практике число конфигураций, используемых в методе MCSCF относительно мало. Тем не менее, использование одних только (квази-вырожденных) конфигураций с наименьшей энергией учитывает статическую корреляцию практически полностью, хотя практически не учитывает динамическую корреляцию.

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

Альтернативой методу MCSCF является теория валентной связи (VB). Теоретически, данный метод порождает столь же великое множество многореферентных методов, как и MCSCF, однако на практике реализована лишь малая их часть. Значительно большее распространение получил «усечённый вариант» метода VB – метод обобщённой валентной связи (GVB).

Данный метод подобен методу HF, но имеет теоретическую «надстройку»: данный метод позволяет попарно объединять внешние молекулярные орбитали так, что внутри каждой пары орбиталей электроны могут свободно распределяться по орбиталям.

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

Существует 2 разновидности метода GVB:

• GVB-PP(n), в котором внешние орбитали объединяются в \(n\) пар, а внутренние описываются с помощью RHF.

• GVB(n), в котором внешние орбитали объединяются в \(n\) пар, а внутренние описываются с помощью ROHF.

Наиболее ярко преимущество GVB выражается в его частном случае – GVB-PP(1), который является полным аналогом двухконфигурационного метода самостогласованного поля (TCSCF), основанного на MCSCF, при этом GVB-PP(1) значительно менее затратный, чем TCSCF.

Метод Кона-Шэма

В отличие от метода Хартри-Фока (HF) для описания всей системы метод Кона-Шэма (KS) использует электронную плотность системы \(\rho\) вместо волновой функции системы \(\Psi\).

Однако, данный метод не отказывается от понятия «орбиталь». В методе Кона-Шэма электронная плотность системы \(\rho\) определяется через орбитали \(\Psi_k\) следующим образом:

\[\rho=\sum^n_{k=1}|\Psi_k|^2\]

При этом спин-орбитали \(\psi_k\) выражаются через базисные функции точно таким же способом, как и в методе Хартри-Фока.

В целом, метод Кона-Шэма (KS) – это лишь одна из реализаций теории функционала плотности (DFT), использующая подходы метода самосогласованного поля (SCF) и понятие спин-орбитали.

Поэтому данный метод также известен как KS-DFT [Kohn-Sham Density Functional Theory]. Существуют и другие методы, например OF-DFT [Orbital-Free Density Functional Theory], в котором отсутстувуют спин-орбитали и понятие волновой функции в целом. Однако метод KS-DFT до сих пор остаётся наиболее проработанной и часто используемой реализацией DFT.

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

Примечание

Функционал – это математический объект, сопоставляющий функции число. Простейщий пример функционала – это значение функции в фиксированной точке.

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

Энергия системы \(E\) в методе Кона-Шэма представляется в виде суммы кинетической энергии \(E_T\), энергии взаимодействия электронов с ядрами \(E_V\), энергии взаимодействия электронов друг с другом \(E_J\) и обменно-корреляционной энергии \(E_{XC}\):

\[E=E_T+E_V+E_J+E_{XC}\]

Для функционалов первых трёх слагаемых (\(E_T\), \(E_V\) и \(E_J\)) известно точное выражение, однако функционал последнего (\(E_{XC}\)), за исключением тривиальных случаев, может быть выражен только приближённо.

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

Стоит отметить, что, ввиду наличия в методе Кона-Шэма спин-орбиталей, существуют его разновидности, аналогичные таковым для метода Хартри-Фока, а именно – RKS, ROKS и UKS.