Моделирование

Ab initio методы расчета

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

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

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

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

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

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

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

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

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

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

Электронное строение системы (молекулы или кристалла) определяется стационарным уравнением Шредингера

\[ <br />
\widehat H\Psi'=E\Psi' .<br />
 \](1)
$ \Psi' $ - волновая функция электронов и ядер, $ \widehat H $ - полный гамильтониан системы.
\[ <br />
\widehat H=\widehat H_e+\widehat H_N+\widehat H_{eN} ,<br />
 \]
где
\[ <br />
\widehat H_e=-\frac12\sum_i\nabla^2(i)+\frac12\sum_{i\ne j}<br />
\frac1{|\vec r_i-\vec r_j|}<br />
 \]
- гамильтониан электронной подсистемы, суммирование проводится по всем электронам;
\[ <br />
\widehat H_N=-\frac12\sum_n\frac1{m_n}\nabla^2(n)<br />
 \]
- гамильтониан ядерной подсистемы, суммирование ведется по всем ядрам;
\[ <br />
\widehat H_{eN}=-\sum_n \left[\sum_i\frac{z_n}{|\vec r_i-\vec<br />
R_n|}+\frac12\sum_{n'\ne n}\frac{z_nz_{n'}}{|\vec R_n-\vec<br />
R_{n'}|}\right]<br />
 \]
- гамильтониан электронно-ядерного взаимодействия; $ z_n $, $ m_n $ и $ \vec R_n $ - заряд, масса и радиус-вектор $ n $-го ядра, $ \vec r_i $ - радиус-вектор $ i $-го электрона,
\[ <br />
\nabla\equiv\vec e_x\frac{\partial}{\partial x}+ \vec e_y<br />
\frac{\partial}{\partial y}+\vec e_z \frac{\partial}{\partial z} ,<br />
 \]
суммирование ведется по всем ядрам и электронам. Здесь и далее используются атомные единицы: $ \hbar=1 $, $ m_e=1 $, $ e=1 $. Волновая функция $ \Psi' $ зависит от спиновых и пространственных координат всех ядер и электронов, и даже для простой системы точно решить уравнение (1) невозможно, поэтому используют ряд
приближений.

Адиабатическое приближение

В адиабатическом приближении, известном еще как приближение Борна-Оппенгеймера, делается следующее предположение. Поскольку масса ядра $ m_n $ много больше массы электрона $ m_e $, средняя скорость
движения ядер много меньше средней скорости движения электронов, то считают, что при достаточно медленном движении ядер электроны всегда успевают за таким движением. Тогда в уравнении (1)
движения электронов и ядер разделяются, и полная волновая функция системы $ \Psi'(q,Q) $ ($ q $ - набор электронных координат, $ Q $ -набор ядерных координат) ищется в виде произведения электронной и ядерной волновых функций

\[ <br />
\Psi'(q,Q)=X(Q)\Psi(q,Q) ,<br />
 \]
где $ \Psi(q,Q) $ - решение уравнения Шредингера для электронов, для фиксированных координат ядер:
\[  \left[\widehat H_e+\widehat<br />
H_{eN}\right]\Psi(q,Q)=E(Q)\Psi(q,Q) .<br />
 \](2)

Одноэлектронное приближение

В многоэлектронной системе все электроны взаимодействуют между собой, движение каждого электрона определяется движением всех остальных. Согласно же одноэлектронному приближению, впервые предложенном Хартри, электроны движутся независимо друг от друга, и движение каждого электрона определяется его координатами в потенциальном поле, создаваемом ядрами и остальными электронами. Тогда движения электронов разделяются, и для каждого электрона может быть введена одноэлектронная функция, являющаяся решением уравнения Шредингера вида (2), где вместо $ |\vec r_i-\vec r_j| $ в $ \widehat H_e $ записан эффективный потенциал $ V_{эфф} $. Волновая функция многоэлектронной системы представляется в виде произведения одноэлектронных функций (орбиталей)

\[ <br />
\Psi=C\prod_{i=1}^n\psi_i .<br />
 \]

Приближение Хартри-Фока

В приближении Хартри не учитывается спин электрона. С учетом принципа Паули многоэлектронная волновая функция должна быть антисимметрична относительно взаимных перестановок любых пар электронов. Этому условию можно удовлетворить, если представить $ \Psi $ в виде определителя Слейтера:

\[ <br />
\small{\Psi=C\vmatrix\psi_1(1)\alpha(1)&\psi_1(1)\beta(1)&\ldots&<br />
\psi_n(1)\alpha(1)&\psi_n(1)\beta(1)\\<br />
\psi_1(2)\alpha(2)&\psi_1(2)\beta(2)&\ldots&<br />
\psi_n(2)\alpha(2)&\psi_n(2)\beta(2)\\<br />
\vdots&\vdots&\ddots&\vdots&\vdots\\<br />
\psi_1(N)\alpha(N)&\psi_1(N)\beta(N)&\ldots&<br />
\psi_n(N)\alpha(N)&\psi_n(N)\beta(N)\\ \endvmatrix,}<br />
 \](3)
где $ N $ - число электронов, $ n=N/2 $ - число орбиталей. Для ортонормированных одноэлектронных волновых функций нормировочный множитель $ C $ равен $ (N!)^{-1/2} $. В общем случае $ \Psi $ должна включать линейную комбинацию определителей типа (3). В случае замкнутой оболочки решение уравнения (2) можно искать в виде одного определителя Слейтера.

Для систем, содержащих неспаренные электроны на однократно заселенных орбиталях $ \psi $, этого сделать не удается. В таких случаях многоэлектронная волновая функция должна быть представлена в виде линейной комбинации определителей Слейтера $ \widetilde\Psi_l $, соответствующих различным вариантам заполнения орбиталей - электронным конфигурациям:

\[ <br />
\Psi=\sum_l\widetilde\Psi_l.<br />
 \]

Для замкнутых оболочек орбитали $ \psi_i(r) $ (ортонормированные) являются решением уравнения Хартри-Фока

\[ <br />
\widehat f_i \psi_i(\vec{r})=\varepsilon_i\psi_i(\vec{r})<br />
 \](4)
с оператором Фока
\[ <br />
\widehat F\equiv \sum_{i=1}^N \equiv \widehat h+\widehat<br />
J-\widehat K=\widehat h+\sum_{l=1}^n \left(2\widehat J_l-\widehat<br />
K_l\right) .<br />
 \]

Собственные значения $ \varepsilon_i $ представляют собой орбитальные энергии. Оператор $ \widehat h $ описывает кинетическую энергию электронов и их взаимодействие с остовом

\[ <br />
\widehat h\equiv-\frac12\sum_1\nabla^2(i)+\sum_i\sum_n\frac{z_n}<br />
{|\vec r_i-\vec R_n|} .<br />
 \](5)

Кулоновские ($ \widehat J_l, \widehat J  $) и обменные ($ \widehat K_l, \widehat K $) операторы определены равенствами

\[ <br />
\widehat J\psi(\vec r)\equiv\left(\int\frac{\rho (\vec r'|\vec<br />
r')}{|\vec r-\vec<br />
r'|}dV'\right)\psi(\vec r),\\<br />
\widehat K\psi(\vec r)\equiv\left(\int\frac{\rho (\vec r|\vec<br />
r')}{|\vec r-\vec r'|}dV'\right)\psi_l(\vec r) ,<br />
 \]

\[ 
\widehat J_l\psi(\vec r)\equiv\left(\int\frac{\psi_l(\vec r')
\psi_l^*(\vec r')}{|\vec r-\vec
r'|}dV'\right)\psi(\vec r),\\
\widehat K_l\psi(\vec r)\equiv\left(\int\frac{\psi(\vec r')
\psi_l^*(\vec r')}{|\vec r-\vec r'|}dV'\right)\psi_l(\vec r) ,
 \]

где

\[ <br />
\rho (\vec r|\vec r')=2\sum_{i=1}^{n}\psi_i (\vec r)\psi_i^* (\vec<br />
r')<br />
 \]
- бесспиновая матрица плотности первого порядка. Кулоновский оператор $ \widehat J_l $ описывает усредненное электростатическое поле, создаваемое всеми электронами кроме $ l $-го. Оператор $ \widehat K_l $ описывает так называемое обменное взаимодействие, не имеющее классического аналога и являющееся следствием учета
принципа Паули. Вводя обозначения
\[ <br />
J_{ij}\equiv\int\!\!\!\int\frac{\psi_i(\vec r)\psi_i^*(\vec r)<br />
\psi_j(\vec r')\psi_j^*(\vec r')}<br />
{|\vec r-\vec r'|}dVdV', \\<br />
K_{ij}\equiv\int\!\!\!\int\frac{\psi_i(\vec r)\psi_j^*(\vec r)<br />
\psi_j(\vec r')\psi_i^*(\vec r')} {|\vec r-\vec r'|}dVdV' ,<br />
 \]
получаем выражение для полной энергии системы
\[ <br />
E=2\sum_{i=1}^n\varepsilon_i-\sum_{j=1}^n\left(2J_{ij}-K_{ij}\right).<br />
 \](6)

Одноэлектронные волновые функции $ \psi_i $ в определителе Слейтера (3), молекулярные орбитали, строятся в виде разложения по базисным функциям $ \phi_{\mu} $:

\[ <br />
\psi_i=\sum_{\mu=1}^K C_{\mu}^i\phi_{\mu}.<br />
 \](7)
Волновые функции $ \psi_i $ находятся вариационным методом путем минимизации полной энергии. Применительно к данной задаче необходимо найти набор коэффициентов $ \{C_{\mu}^i\} $. Для этого необходимо решить систему уравнений Рутаана-Холла (сам метод получил название Хартри-Фока-Рутаана):
\[ <br />
\sum_{\nu=1}^K (F_{\mu\nu}-\varepsilon_iS_{\mu\nu})\cdot<br />
C_{\mu}^i=0, \quad (\mu = 1, 2, ..., K),<br />
 \](8)
где $ K $ - размер базиса, $ S_{\mu\nu} $ - интегралы перекрывания:
\[ <br />
S_{\mu\nu}=\int\phi_{\mu}^*(1)\phi_{\nu}(1)dV_2,<br />
 \](9)
$ F_{\mu\nu} $ - матричный элемент одноэлектронного оператора Фока
$ \widehat F $:
\[ <br />
F_{\mu\nu}=h_{\mu\nu}+\sum_{\lambda=1}^K\sum_{\sigma=1}^K<br />
P_{\lambda\sigma}[(\mu\nu|\lambda\sigma)<br />
-\frac12(\mu\lambda|\nu\sigma)],<br />
 \](10)
включающий матричные элементы опрератора $ \widehat h $ (5), кулоновский и обменный члены. $ P_{\mu\nu} $ - элементы матрицы плотности:
\[ <br />
 P_{\mu\nu}=2\sum_{i=1}^{occup.}C_{\mu }^{i*} \cdot C_{\nu}^i,<br />
 \](11)
суммирование идет по всем занятым орбиталям, $ (\mu\nu|\lambda\sigma) $ и $ (\mu\lambda|\nu\sigma) $ -
двухцентровые интегралы:
\[ <br />
 (\mu\nu|\lambda\sigma)=\int\int\phi_{\mu}^*(1)\phi_{\nu}(1)\frac12<br />
 \phi_{\lambda}^*(2)\phi_{\sigma}(2)dV_1dV_2,<br />
 \]

Отметим, что для каждого $ i $ мы имеем систему из $ K $ уравнений (8), содержащих $ K $ неизвестных $ \{C_{\mu}^i\} $. Однозначность решения системы уравнений (8) обеспечивает условие нормировки одноэлектронных волновых функций $ \psi_i $:

\[ <br />
 \int\psi_i^*\psi_idV=1=\sum_{\mu=1}^K\sum_{\nu=1}^K<br />
C_{\mu}^{i*}S_{\mu\nu} C_{\nu}^{i}.<br />
 \]

Можно заметить, что для нахождения коэффициентов $ \{C_{\mu}^i\} $ необходимо построить матрицу Фока $ F $, для чего необходимо знать коэффициенты $ \{C_{\mu}^i\} $. Для решения задачи на первом этапе задают значения коэффициентов $ \{C_{\mu}^i\} $, иными словами выбирают начальное приближение. После чего Решают систему
уравнений (8) и находят новые значения для $ \{C_{\mu}^i\} $. Если разница между новыми и старыми значениями коэффициентов меньше некоторого заранее заданного параметра $ \delta $, то процедура повторяется до тех пор, пока не будет выполнено условие:

\[ <br />
 |\{C_{\mu}^i\}^{n} - \{C_{\mu}^i\}^{n-1}| < \delta.<br />
 \]

Такая процедура называется процедурой самосогласованного поля (self-consistent field, SCF). В ходе процедуры самосогласования полная энергия системы понижается, таким образом достижение сходимости по энергии выступает в роли дополнительного и надежного критерия завершения SCF процедуры.

A posteriori учет корреляционных эффектов

Недостатком метода Хартри-Фока является пренебрежение корреляционными эффектами (приближение Хартри). На самом деле движение электронов не является независимым друг от друга, а скореллировано. Энергию корреляции принято определять как разницу между "точной" энергией системы и энергией, полученной в методе Хартри-Фока $ E_{\textrm{HF}} $:

$ <br />
E_{\textrm{corr}} = E - E_{\textrm{HF}}.<br />
 $

Поскольку энергия Хартри-Фока всегда выше "точной" энергии, энергия корреляции всегда отрицательна:

$ E_{\textrm{corr}} < 0.<br />
 $

Для учета эффектов корелляции разработан ряд a posteriori методов. Рассмотрим наиболее популярные из них.

Метод Меллера-Плессе

Теория возмущений Меллера-Плессе (Møller-Plesset perturbation theory) одна из наиболее популярных a posteriory коррекций метода Хартри-Фока, позволяющая учесть эффекты корреляции электронов. Основная идея метода была опубликована еще в 1934 году.

Рассмотрим невозмущенный Гамильтониан $ \widehat H_0 $ к которому добавлено небольшое возмущение $ \widehat V $:

\[ <br />
\widehat H = \widehat H_0 + \lambda\widehat V ,<br />
 \](1)
где $ \lambda $ есть некий параметр (безразмерная величина). В теории Меллера-Плессе нулевого порядка волновые функции являются собственными функциями оператора Фока, который выступает в роли невозмущенного оператора. Возмущение в данном случае есть корреляционный потенциал.

Представим возмущенную волновую функцию и соответствующую ей энергию в виде степенного ряда по $ \lambda $:

\[ <br />
\Psi = \lim_{n \to \infty}\sum_{i=0}^n\lambda^i\Psi^{(i)},<br />
 \](2)
\[ <br />
E =\lim_{n \to \infty}\sum_{i=0}^n\lambda^i E^{(i)}.<br />
 \](3)

Подставим теперь (2) и (3) в стационарное уравнение Шредингера ($ n \to \infty $):

\[ <br />
\left(\widehat H_0 + \lambda\widehat V\right)<br />
\left(\sum_{i=0}^n\lambda^i\Psi^{(i)}\right) =<br />
\left(\sum_{i=0}^n\lambda^i<br />
E^{(i)}\right)\left(\sum_{i=0}^n\lambda^i\Psi^{(i)}\right),<br />
 \](4)

Приравнивание множителей при $ \lambda^k $ дает $ k $-тый порядок теории возмущений ($ k = 1, 2, ..., n $):

\[ <br />
(\widehat H_0 - E^{(0)}) \Psi^{(0)}= 0<br />
 \](5)
\[ <br />
(\widehat H_0 -E^{(0)})\Psi^{(1)}= (E^{(1)}-V) \Psi^{(0)}<br />
 \](6)
\[ <br />
(\widehat H_0 - E^{(0)})\Psi^{(2)}= (E^{(1)}-V)\Psi^{(1)}<br />
+E^{(2)}\Psi^{(0)}<br />
 \](7)
и т.д.

Возмущающий, в данном случае корреляционный, потенциал $ V $ есть:

\[ <br />
\widehat V \equiv \widehat H - \widehat F = \widehat H -<br />
\sum_k\widehat f_k,<br />
 \](8)

где $ \Psi_0 $ - нормализованный определитель Слейтера, который есть собственная функция $ \widehat H_0 $.

Теперь будем поочередно рассматривать уравнения (5), (6), (7), что даст нам выражения для поправок к энергии $ E^{(i)} $.

В нулевом порядке теории возмущений имеем:

\[ <br />
E^{(0)} = \sum_k \varepsilon_k,<br />
 \](9)
что есть сумма энергий орбиталей.

Поправка первого порядка определяется выражением

\[ <br />
E^{(1)} = \langle \Psi^{(0)}|V|\Psi^{(0)}\rangle.<br />
 \](10)

Если мы остановимся на учете только первой поправки, то получим

\[ <br />
E= E^{(0)} + E^{(1)} = \langle \Psi^{(0)}|H_0+V|\Psi^{(0)}\rangle,<br />
 \](11)
что есть энергия метода Хартри-Фока.

Чтобы вычислить поправку $ E^{(2)} $, необходимо найти сначала волновую функцию $ \Psi^{(1)} $:

\[ <br />
\Psi^{(1)}=\sum_t\frac{|\langle<br />
\Psi^{(0)}|V|\Psi_t\rangle|}{E_t-E^{(0)}}, \ \ \ \ (t>0),<br />
 \](12)
после чего можно записать вторую поправку к энергии как:
\[ <br />
E^{(2)}=-\sum_t\frac{|\langle<br />
\Psi^{(0)}|V|\Psi_t\rangle|^2}{E_t-E^{(0)}} \cdot \Psi_t, \ \ \ \<br />
(t>0).<br />
 \](13)

Отметим, что $ E^{(2)}<0 $, поэтому MP2 метод всегда приводит к понижению энергии. Однако, поскольку метод не является вариационным, то MP2 может занижать истинное значение энергии и последующие поправки (MP3, MP4, ...) будут давать более высокие значения энергии, по сравнению с MP2.

Метод конфигурационного взаимодействия

Другим методом, позволяющим учесть электронные корреляции, является метод конфигурационного взаимодействия(сonfiguration interaction, CI). С точки зрения математики конфигурация описывает линейную комбинацию определителей Слейтера. При описании заселенности орбиталей (например, $ 1s^22s^2p^1 $) взаимодействие означает смешивание различных электронных конфигураций (состояний). Однако из-за большого компьютерного времени, требуемого для выполнения CI расчетов, обычно этот метод применяется для небольших систем.

Итак, напомним вариационный принцип:

\[ <br />
E_0 \equiv \underbrace{min}_{\Psi} \langle\Psi|\widehat<br />
H|\Psi\rangle. <br />
 \](1)
(предполагается, что функции нормированы). Мы можем получить "точную" энергию основного состояния и "точную" волновую функцию, учитывая корреляционные эффекты, минимизируя энергию по всем норированным антисимметричным волновым функциям. Однако это задача очень сложна. В самом простом случае мы можем рассмотреть один определитель Слейтера, что приводит к уравнениям Хартри-Фока. Однако, в случае, когда различные определители Слейтера имеют близкие энергии (например для молекулы озона), рассмотрение одного определителя Слейтера является плохим приближением. В этом случае необходимо рассматривать вариационную волновую функцию, которая является линейной комбинацией функций, соответствующих различным конфигурациям:
\[ <br />
\Psi = \sum_k c_k\Psi_k.<br />
 \](2)

Напомним, что для заданного оператора в эрмитовом пространстве набор собственных функций операторя образует полный ортонормированный набор. Это означает, что если мы рассматриваем определители Слейтера, то разложение (2) будет точным. Это аналогично разложению в теории возмущений. В методе CI разложение идет по всем возможным определителям Слейтера, составленным из хартри-фоковских орбиталей.

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

Для удобства сумму (2) можно переписать, сгруппировав определители Слейтера с одним "возбужденным" электроном $ \Psi_i^a $ (электрон перешел с $ i $-той занятой орбитали на $ a $-тую незанятую), с двумя "возбужденными" электронами $ \Psi_{ij}^{ab} $ и т.д. Тогда "точная" волновая функция будет иметь вид

\[ <br />
\Psi = c_0\Psi_0 + \sum_{i=1}^N\sum_{a=N+1}^K c_i^a\Psi_i^a +<br />
\sum_{i>j=1}^N\sum_{a>b=N+1}^K c_{ij}^{ab}\Psi_{ij}^{ab} +<br />
 \](3)
$$<br />
+ \sum_{i>j>k=1}^N\sum_{a>b>c=N+1}^K c_{ijk}^{abc}\Psi_{ijk}^{abc} +<br />
\sum_{i>j>k>l=1}^N\sum_{a>b>c>d=N+1}^K<br />
c_{ijkl}^{abcd}\Psi_{ijkl}^{abcd} + ...<br />
$$

Заметим, что $ \Psi_{ij}^{ab}=-\Psi_{ji}^{ab}=-\Psi_{ij}^{ba}=\Psi_{ji}^{ba} $.

На практике невозможно рассчитать все незанятые орбитали. Также совершенно очевидно, что число членов в разложении очень быстро возрастает при учете все более и более высоких уровней возбуждения.

Для удобства записи введем вектора $ \textbf{\textit{i}} $ и $ \textbf{\textit{a}} $, чьими компонентами являются занятые ($ i_1,i_2, i_3, ...  $) и незанятые ($ a_1, a_2, a_3, ...  $) орбитали, соответственно:

\[ 
\Psi =
c_{[1,2,3...N]}^{[1,2,3...N]}\Psi_{[1,2,3...N]}^{[1,2,3...N]} +
\sum_{i_1=1}^N\sum_{a=N+1}^K
c_{[i_1,i_2...i_N]}^{[a_1,i_2...i_N]}\Psi_{[i_1,i_2...i_N]}^{[a_1,i_2...i_N]}+
 \](4)
$$
+\sum_{i_1>i_2=1}^N\sum_{a_1>a_2=N+1}^K
c_{[i_1,i_2,i_3...i_N]}^{[a_1,a_2,i_3...i_N]}
\Psi_{[i_1,i_2,i_3...i_N]}^{[a_1,a_2,i_3...i_N]} +
\sum_{i_1>i_2>i_3=1}^N\sum_{a_1>a_2>a_3=N+1}^K
c_{[i_1,i_2,i_3,i_4...i_N]}^{[a_1,a_2,a_3,i_4...i_N]}
\Psi_{[i_1,i_2,i_3,i_4...i_N]}^{[a_1,a_2,a_3,i_4...i_N]} +
$$
$$
+\sum_{i_1>i_2>i_3>i_4=1}^N\sum_{a_1>a_2>a_3>i_4=N+1}^K
c_{[i_1,i_2,i_3,i_4,i_5...i_N]}^{[a_1,a_2,a_3,a_4,i_5...i_N]}
\Psi_{[i_1,i_2,i_3,i_4,i_5...i_N]}^{[a_1,a_2,a_3,a_4,i_5...i_N]} + ...
=\sum_{\textbf{\textit{i}},\textbf{\textit{a}}}c_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}
\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}
$$

Для расчета CI волновой функции начнем с уравнения Шредингера

\[ <br />
\widehat H \Psi = E\Psi,<br />
 \](6)
затем домножим обе стороны равенства слева на определитель Слейтера $ \Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}} $:
\[ <br />
\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}\widehat H \Psi =<br />
\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}E\Psi<br />
 \](6)
и, подставив (4) в (6), имеем
\[ <br />
\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}\widehat H<br />
\sum_{\textbf{\textit{j}},\textbf{\textit{b}}}c_{\textbf{\textit{j}}}^{\textbf{\textit{b}}}<br />
\Psi_{\textbf{\textit{j}}}^{\textbf{\textit{b}}} =<br />
E\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}\sum_{\textbf{\textit{j}},\textbf{\textit{b}}}c_{\textbf{\textit{j}}}^{\textbf{\textit{b}}}<br />
\Psi_{\textbf{\textit{j}}}^{\textbf{\textit{b}}}.<br />
 \](7)
После интегрирования имеем:
\[ <br />
\sum_{\textbf{\textit{j}},\textbf{\textit{b}}}\langle\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}<br />
|\widehat H |\Psi_{\textbf{\textit{j}}}^{\textbf{\textit{b}}}<br />
\rangle c_{\textbf{\textit{j}}}^{\textbf{\textit{b}}}=<br />
E\sum_{\textbf{\textit{j}},\textbf{\textit{b}}}\langle\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}<br />
|\Psi_{\textbf{\textit{j}}}^{\textbf{\textit{b}}} \rangle<br />
c_{\textbf{\textit{j}}}^{\textbf{\textit{b}}}  \nonumber\\<br />
=<br />
E\sum_{\textbf{\textit{j}},\textbf{\textit{b}}}\delta_{\textbf{\textit{ij}}}\delta_{\textbf{\textit{ab}}}<br />
c_{\textbf{\textit{j}}}^{\textbf{\textit{b}}} \\<br />
= E c_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}.\nonumber<br />
 \](8)
Это есть ничто иное как задача на нахождение собственных значений матрицы гамильтониана $ \widehat H $ с матричными элементами
\[ <br />
H_{\textbf{\textit{i,j}}}^{\textbf{\textit{a,b}}}\equiv \langle<br />
\Psi_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}|\widehat<br />
H|\Psi_{\textbf{\textit{j}}}^{\textbf{\textit{b}}} \rangle.<br />
 \](9)
Нахождение наименьшего собственного значения из
\[ <br />
\sum_{\textbf{\textit{i}},{\textbf{\textit{a}}}}<br />
H_{\textbf{\textit{i,j}}}^{\textbf{\textit{a,b}}}<br />
c_{\textbf{\textit{j}}}^{\textbf{\textit{b}}} = E<br />
c_{\textbf{\textit{i}}}^{\textbf{\textit{a}}}<br />
 \](10)
равносильно минимизации энергии по всем волновым функциям вида (4).

Естественно, что на практике разложение (4) не рассматривают выше какого-либо порядка. Если оборвать разложение после нулевого порядка, то имеем метод Хартри-Фока. Если после первого, то имеем Конфигаруционное Взаимодействие с Одноэлектронными Возбуждениями (Configuration-Interaction with Single excitations, CIS), если после второго порядка, то Конфигаруционное Взаимодействие с Одно и Двухэлектронными возбуждениями (Configuration Interaction with Single and Double excitations, CISD) и т.д.: (CISDT - третьего порядка, CISDTQ - четвертого порядка) ... Если разложение не обрывать и учитывать все $ N $-электронные возбуждения, то имеем метод Полного Конфигурационного Взаимодействия Full-Configuration-Interaction, FCI), однако этот метод требует огромных вычислительных ресурсов, которые можно оценить с помощью биноминального коэффициента

\[ <br />
\begin{pmatrix}<br />
  K\\<br />
  N\\<br />
\end{pmatrix}\equiv \frac{(K)!}{(K-N)!N!},<br />
 \](11)
где $ K $ - полное число орбиталей хартри-фоковских орбиталей (занятых и свободных), $ N $ - число электронов в системе. Для больших $ K $ FCI расчеты достаточно точны, и точность других методов расчета может быть оценена в рамках критерия "насколько близки полученные результаты к результатам FCI".

Метод связанного кластера

Другим методом, основанным на методе Хартри-Фока и позволяющим учесть электронные корреляции, является метод Связанного кластера (Coupled Cluster, CC).

Изначально метод был разработан Fritz Coester и Hermann Kummel еще в 1959-х годах для исследования явлений, связанных с ядерной физикой. Однако наиболее широкое распространение получил в 1960-х годах, после того, как был переформулирован Jiri Cizek и Josef Paldus для расчета энергии корреляции атомов и молекул.

В методе CC решение уравнения Шредингера $ \Psi $ представляется в следующем виде:

\[ <br />
\Psi=e^{\widehat T}\Psi_0,<br />
 \](1)
где $ \Psi_0 $ - определитель Слейтера, $ \widehat T $ - оператор возбуждения, который, действуя на $ \Psi_0 $ производит линейную комбинацию возбужденных детерминантов Слейтера.

Кластерный оператор $ \widehat T $ может быть представлен в виде ряда

\[ <br />
\widehat T = \widehat T_1 + \widehat T_2 +\widehat T_3 + ... ,<br />
 \](2)
где $ \widehat T_1 $ - оператор одноэлектронного возбуждения, $ \widehat T_2 $ - оператор двухэлектронного возбуждения и т.д. имеют вид
\[ <br />
\widehat T_1 = \sum_i\sum_a t_i^a \widehat a_i \widehat<br />
a_a^{\dagger},\nonumber\\<br />
\widehat T_2 = \frac14\sum_{i,j}\sum_{a,b} t_{ij}^{ab} \widehat<br />
a_i \widehat a_j\widehat a_a^{\dagger}\widehat a_b^{\dagger}<br />
 \](3)
и т.д.

В формуле (3) $ \widehat a^{\dagger} $ и $ \widehat a $ обозначают опрераторы рождения и уничтожения, соответственно. Индексы $ i,j $ и $ a,b $ как и ранее нумеруют занятые и незанятые орбитали.

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

Для нахождения $ \Psi $ необходимо найти коэффициенты разложения $ t_i^a $ и $ t_{ij}^{ab} $. Учитывая вид оператора $ \widehat T $, разложим экспоненциальный оператор $ e^{\widehat T} $ в ряд Тейлора:

\[ <br />
e^{\widehat T} = 1+ \widehat T +\frac{\widehat T^2}{2!} + ... =<br />
1+\widehat T_1+\widehat T_2 +\frac{\widehat T_1^2}{2!} + \widehat<br />
T_1\widehat T_2 + \frac{\widehat T_2^2}{2!} +...<br />
 \](4)

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

Теория функционала плотности

Альтернативой методу Хартри-Фока является метод, основывающийся на теории функционала плотности Density functional theory - DFT). Основным преимуществом данного метода является то, что эффекты корреляции удается учесть сразу, не пребегая к a posteriori коррекции, что позволяет существенно сэкономить время
счета.

В основе метода лежит предположение о том, что полная энергия есть функция плотности $ \rho(\vec{r}) $. Это означает, что нет необходимости знать сложную многоэлектронную волновую функцию.

Согласно теореме Хоэнберга-Кона, плотность $ \rho(\vec{r}) $ для основного состояния многоэлектроной системы в некотором внешнем потенциале однозначно определяет этот потенциал. А значит, $ \rho(\vec{r}) $ неявно определяет все свойства системы, которые можно получить решая уравнение Шредингера. Исходя из этих рассуждений функционал полной энергии $ E_{tot} $ можно представить в виде

\[ 
E_{tot}[\rho(\vec{r})] \equiv T_0[\rho(\vec{r})]+\int V_0(\vec{r})\rho(\vec{r})d\vec{r}
+\frac{1}{2}\int\frac{\rho(\vec{r})\rho(\vec{r'})}{|\vec{r}-\vec{r'}|}d\vec{r}d\vec{r'}
+E_{xc}[\rho(\vec{r})], \ \quad
 \](7)

где первое слагаемое -- кинетическая энергия невзаимодействующих электронов, второе -- взаимодействие между электронами и ядрами, задаваемое потенциалом $ V_0(\vec{r}) $ и электронной плотностью $ \rho(\vec{r}) $, третье -- кулоновское взаимодействие между электронами. Наконец последний член содержит обменно-корреляционную энергию $ E_{xc}[\rho(\vec{r})] $, которая определяет все остальные вклады в межэлектронное взаимодействие.

Записав функционал полной энергии $ E_{tot} $ в виде (7) можно применить к нему вариационный принцип и получить одноэлектронные уравнения Кона-Шэма:

$ <br />
\{-\frac{1}{2}\nabla^2+V_{eff}(\vec{r})\}\Phi_i(\vec{r})=\varepsilon_i\Phi_i(\vec{r}),<br />
 $

где
$ <br />
   V_{eff}(\vec{r}) = V_0(\vec{r})+V_C[\rho(\vec{r})]+\mu_{xc}[\rho(\vec{r})]<br />
 $

есть эффективный потенциал в котором находятся электроны, а $ \mu_{xc} $ есть производная функционала $ E_{xc} $ по плотности:

\[ 
  \mu_{xc}[\rho(\vec{r})] = \frac{\delta E_{xc}[\rho(\vec{r})]}{\delta \rho(\vec{r})} \equiv
   \frac{\delta\{\int \rho(\vec{r'})\varepsilon_{xc}[\rho(\vec{r'})]d\vec{r'}\}}{\delta \rho(\vec{r})}.
 \]

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

$ <br />
\rho(\vec{r})=\sum_{i\ occupied}|\psi_i(\vec{r})|^2.<br />
 $

Отметим, что в отличие от метода Хартри-Фока собственные функции $ \psi_i(\vec{r}) $ и собственные значения $ \varepsilon_i $ уравнений Кона-Шэма не имеют строгого физического смысла.

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

Приближение локальной плотности

В рамках приближение локальной плотности (LDA) изначально неоднородная система трактуется как локально однородная (аналог электронного газа) с обменно-корреляционной энергией $ \varepsilon_{xc}[\rho(\vec{r})] $, которая определяется выражением

$ <br />
E_{xc}^{LDA}[\rho(\vec{r})] = \int<br />
\rho(\vec{r})\varepsilon_{xc}[\rho(\vec{r})] d\vec{r}.<br />
 $

Для магнитных систем, где допускается спиновая поляризация, вводится приближение локальной спиновой плотности (LSDA), и обменно-корреляционная энергия $ E_{xc} $ есть функционал локальных электронно-спиновых плотностей $ \rho_{\uparrow} $ и $ \rho_{\downarrow} $:

$ <br />
E_{xc}^{LSDA}[\rho_{\uparrow}(\vec{r}),<br />
\rho_{\downarrow}(\vec{r})] = \int<br />
\rho(\vec{r})\varepsilon_{xc}[\rho_{\uparrow}(\vec{r}),<br />
\rho_{\downarrow}(\vec{r})] d\vec{r},<br />
 $

Как правило $ E_{xc}[\rho(\vec{r})] $ представляют в виде суммы обменного и корреляционного вкладов:

$ <br />
E_{xc}[\rho(\vec{r})] = E_{x}[\rho(\vec{r})]+E_{c}[\rho(\vec{r})].<br />
 $

Наиболее известные обменные LDA функционалы

К одним из наиболее ранних версий LDA-DFT относится так называемый метод $ X_{\alpha} $, в котором обменный функционал имеет вид:

$ <br />
E_{x}^{X_{\alpha}}[\rho] =<br />
-\frac{9\alpha}{4}\left(\frac{3\rho}{\pi}\right)^{\frac13}.<br />
 $
Значение подгоночного параметра $ \alpha $ лежит в диапазоне от $ \frac23 $ до 1. Функционал с $ \alpha=\frac23 $ известен как функционал Дирака, а с $ \alpha=1 $ - как функционал Слейтера.

Наиболее известные корреляционные LDA функционалы

  • Vosko-Wilk-Nusair
  • Perdew-Zunger
  • Cole-Perdew
  • Perdew-Wang

Приближение обобщенного градиента

Несмотря на свою простоту метод LDA часто дает весьма удовлетворительное согласие с экспериментальными результатами. Однако в ряде случаев его оказывается недостаточно. Естественным шагом является градиентная коррекция функционала (приближение обобщенного градиента - GGA), где $ E_{xc} $ рассматривается не только как функция спиновых плотностей, но и их градиентов:

\[ 
E_{xc}^{GGA}[\rho_{\uparrow}(\vec{r}), \rho_{\downarrow}(\vec{r})]
= \int f[\rho_{\uparrow}(\vec{r}), \rho_{\downarrow}(\vec{r}),
\nabla\rho_{\uparrow}(\vec{r}), \nabla\rho_{\downarrow}(\vec{r})]
d\vec{r},
 \](5)

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

В настоящее время большой популярностью пользуются гибридные функционалы, где включен "точный" Хартри-Фоковский обменный член:

\[ <br />
E^{h}_{xc}=aE^{HF}_{xc} + bE^{DFT}_{xc},<br />
 \](6)
где $ a $ и $ b $ - параметры. Так, например, широкой популярностью пользуется гибридный функционал B3LYP.

Расчеты методом DFT обычно выполняются по тем же схемам, что и в методе Хартри-Фока. Различие заключается в расчете обменно-корреляционной энергии, которая не может выть вычислена аналитически и рассчитывается численно.

Среди способов расчета $ E_{xc} $ выделяют сеточные и бессеточные методы. В сеточных расчетах используют сетку точек в пространстве, равномерную (однинаковое количество точек на каждой сфере), либо усеченную (наибольшая плотность точек приходится на те области, в которых свойства системы меняются наиболее быстро)

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

Расчет параметров спектров магнитного резонанса

Расчет g-фактора

Расчет сверхтонких полей

Расчет тензора ГЭП

Расчет химических сдвигов

Молекулярная динамика