怎样写一个量子化学计算程序

开新坑，怎样写一个量子化学程序，慢慢更新，最迟六月份更新完。

这篇文章是基于我的草稿由机器生成，排版、语言和符号可能还有大量错误，我还没有修正。

\(
\def \bra#1{{\langle #1 \rvert}}
\def \ket#1{{\lvert #1 \rangle}}
\def \braket#1#2{{\langle #1 | #2 \rangle}}
\def \expect#1{{\langle #1 \rangle}}
\)

1. 量子化学计算引论

首先，如果用一句话概括量子化学计算，那就是我们在这里要讲到的方法，或者说当代的很多量子化学方法，多数属于变分和微扰的方法。变分方法本质或者说核心问题，在于寻找合适的变分域和线性化方法。

在计算方法上，有迭代的方法和非迭代方法之分，例如， Coupled Cluster 方法本质是迭代的方法，而微扰理论(perturbation theory)和组态相互作用(Configuration Interaction)方法本质上是非迭代的方法，而实际求解CI问题的时候则是用迭代的方法。

\[\text{===Born$–$Oppenheimer approximation===}\]

首先，我们来复习波恩-奥本海默近似的概念。

系统的总哈密顿算符可以被写为如下的形式

\[\hat{H}_{\text{tot}}\text{ = }\hat{T}_e\text{ + }\hat{T}_n\text{+ }\hat{V}_{\text{ee}}\text{ + }\hat{V}_{\text{ne}}\text{ + }\hat{V}_{\text{nn}}\]

系统的薛定谔方程是

\[\hat{H}_{\text{tot}}\left|\psi _{\text{tot}}\text{$\rangle $=}E_{\text{tot}}\right|\psi _{\text{tot}}\rangle\]

解在坐标表象下为

\[\text{$\langle $r,R$|$}\psi _{\text{tot}}\text{$\rangle $=($\langle $r$|\otimes \langle $R$|$)$|$}\psi _{\text{tot}}\text{$\rangle $
}\]

注意：Subscript[Overscript[T, ^], e], Subscript[Overscript[T, ^], n] 不依赖 r, R, 而 Subscript[Overscript[V, ^], ee] 仅依赖 r，Subscript[Overscript[V, ^], nn] 仅依赖 R，Subscript[Overscript[V, ^], ne] 依赖 r，R

定义

\[\hat{H}_e\text{ = }\hat{T}_e\text{ + }\hat{V}_{\text{ee}}\text{ + }\hat{V}_{\text{ne}}\text{ + }\hat{V}_{\text{nn}}\]

\[\text{$\langle $r,R$|$}\hat{H}_{\text{mp}}\text{ = -}\frac{1}{2M_{\text{tot}}}\left(\sum _{i=1}^{N_{\text{elec}}} \nabla _i )^2\right.\]

则有

\[\hat{H}_{\text{tot}}\text{ = }\hat{H}_e\text{ + }\hat{H}_{\text{mp}}\text{ + }\hat{T}_n\]

Subscript[Overscript[H, ^], mp]被称为质量极化 (mass polarization)

由电子薛定谔方程

\[\hat{H}_e\text{(R)$|$}\psi _i\text{(R)$\rangle $ = }E_i\text{(R)$|$}\psi _i\text{(R)$\rangle $}\]

可解得 {\[LeftAngleBracket]r\[VerticalSeparator]Subscript[\[Psi], i](R)\[RightAngleBracket]}，由 Subscript[Overscript[H, ^], e](R) 算符是 Hermitian Operator, {\[LeftAngleBracket]r\[VerticalSeparator]Subscript[\[Psi], i](R)\[RightAngleBracket]} 构成一组完备正交基，从而 Subscript[Overscript[H, ^], tot] 的 Eigen Wavefunction 可以展开在这组完备基上

\[\text{$\langle $r,R$|$}\psi _{\text{tot}}\text{$\rangle $ = ($\langle $R$|\otimes \langle $r$|$)(}\sum _{i=1}^{\infty } |\psi _i\text{(R)$\rangle \langle $}\psi _i\text{(R)$|$)$|$}\psi _{\text{tot}}\text{$\rangle $ = }\sum _{i=1}^{\infty } \text{$\langle $r$|$}\psi _i\text{(R)$\rangle $($\langle $R$|\langle $}\psi _i\text{(R)$|$}\psi _{\text{tot}}\text{$\rangle $)}\]

定义

\[\left.\text{$\langle $R$|\langle $}\psi _i\text{(R)$|$}\psi _{\text{tot}}\text{$\rangle $ = $\langle $R$|$}\psi _{\text{ni}}\right\rangle\]

则

\[\left.\text{$\langle $r,R$|$}\psi _{\text{tot}}\text{$\rangle $ = }\sum _{i=1}^{\infty } \text{$\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\right\rangle\]

现在，我们来写出哈密顿算符的作用

\[\text{$\langle $r,R$|$}\hat{T}_n\text{ = }\sum _{a=1}^{N_{\text{nuc}}} -\frac{1}{2M_a}\nabla _a{}^2\text{ = }\nabla _n{}^2\]

\[\text{$\langle $r,R$|$}\hat{H}_{\text{tot}}|\psi _{\text{tot}}\text{$\rangle $ = $\langle $r,R$|$(}\hat{T}_n\text{ + }\hat{H}_e\text{ + }\hat{H}_{\text{mp}}\text{)$|$}\psi _{\text{tot}}\text{$\rangle $ }\]

\[\left.\text{= (}\nabla _n{}^2\text{ + $\langle $r,R$|$}\hat{H}_e\text{ + $\langle $r,R$|$}\hat{H}_{\text{mp}}\right)\sum _{i=1}^{\infty } \text{$\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ }\]

\[\text{= }\sum _{i=1}^{\infty } \left(\nabla _n{}^2\text{($\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $) + $\langle $r,R$|$}\hat{H}_e\text{($\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $) + $\langle $r,R$|$}\hat{H}_{\text{mp}}\text{($\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $))}\right.\]

\[\text{= }\sum _{i=1}^{\infty } \left(\nabla _n\text{($\langle $R$|$}\psi _{\text{ni}}\right\rangle \nabla _n\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $ + $\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $) + }E_i\text{(R)$\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{ni}}\text{$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $) }\]

\[\left.\text{= }\sum _{i=1}^{\infty } \text{($\langle $R$|$}\psi _{\text{ni}}\right\rangle \nabla _n{}^2\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $ + 2}\nabla _n\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n{}^2\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + }E_i\text{(R)$\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{ni}}\text{$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $)}\]

\[\left.\text{= }E_{\text{tot}} \sum _{i=1}^{\infty } \text{ $\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\right\rangle\]

上式两边左乘 \[LeftAngleBracket]Subscript[\[Psi], j](R)\[VerticalSeparator]r\[RightAngleBracket]，得

\[\left\langle \psi _j\text{(R)$|$r$\rangle $}\sum _{i=1}^{\infty } \text{($\langle $R$|$}\psi _{\text{ni}}\right\rangle \nabla _n{}^2\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $ + 2}\nabla _n\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n{}^2\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + }E_i\text{(R)$\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{ni}}\text{$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $) = $\langle $}\psi _j\text{(R)$|$r$\rangle $ }E_{\text{tot}} \sum _{i=1}^{\infty } \text{ $\langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\rangle\]

\[\left.\sum _{i=1}^{\infty } \text{($\langle $R$|$}\psi _{\text{ni}}\rangle \langle \psi _j\text{(R)$|$r$\rangle $}\nabla _n{}^2\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $ + 2$\langle $}\psi _j\text{(R)$|$r$\rangle $}\nabla _n\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $}\psi _j\text{(R)$|$r$\rangle \langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n{}^2\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + }E_i\text{(R)$\langle $}\psi _j\text{(R)$|$r$\rangle \langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{ni}}\rangle \langle \psi _j\text{(R)$|$r$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $) = }E_{\text{tot}} \sum _{i=1}^{\infty } \text{ $\langle $}\psi _j\text{(R)$|$r$\rangle \langle $r$|$}\psi _i\text{(R)$\rangle \langle $R$|$}\psi _{\text{ni}}\right\rangle\]

由正交关系

\[\left\langle \psi _j\text{(R)$|$}\psi _i\text{(R)$\rangle $ = }\delta _{\text{ij}}\right.\]

得

\[\left.\nabla _n{}^2\text{$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + }E_j\text{(R)$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + }\sum _{i=1}^{\infty } \text{($\langle $R$|$}\psi _{\text{ni}}\rangle \langle \psi _j\text{(R)$|$r$\rangle $}\nabla _n{}^2\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $ + 2$\langle $}\psi _j\text{(R)$|$r$\rangle $}\nabla _n\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $}\nabla _n\text{$\langle $R$|$}\psi _{\text{ni}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{ni}}\rangle \langle \psi _j\text{(R)$|$r$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _i\text{(R)$\rangle $) = }E_{\text{tot}}\text{$\langle $R$|$}\psi _{\text{nj}}\right\rangle\]

上式一共有5项，前2项中已经没有电子波函数，在求和号中的3项里，第1项是2阶非绝热耦合元(second-order non-adiabatic coupling elements)，第2项是1阶非绝热耦合元，而最后一项则是质量极化。

所谓的绝热近似(adiabatic approximation)，就是指将总波函数限制在一个 electronic surface 上。因此，上式中所有非对角元都被忽略，我们得到

\[\left.\nabla _n{}^2\text{$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + }E_j\text{(R)$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{nj}}\rangle \langle \psi _j\text{(R)$|$r$\rangle $}\nabla _n{}^2\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $ + 2$\langle $}\psi _j\text{(R)$|$r$\rangle $}\nabla _n\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $}\nabla _n\text{$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{nj}}\rangle \langle \psi _j\text{(R)$|$r$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $ = }E_{\text{tot}}\text{$\langle $R$|$}\psi _{\text{nj}}\right\rangle\]

注意到

\[\left\langle \psi _j\text{(R)$|$r$\rangle $}\nabla _n\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $ = 0}\right.\]

因此

\[\left.\nabla _n{}^2\text{$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + }E_j\text{(R)$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ + $\langle $R$|$}\psi _{\text{nj}}\rangle \langle \psi _j\text{(R)$|$r$\rangle $}\nabla _n{}^2\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $ + $\langle $R$|$}\psi _{\text{nj}}\rangle \langle \psi _j\text{(R)$|$r$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $ = }E_{\text{tot}}\text{$\langle $R$|$}\psi _{\text{nj}}\right\rangle\]

或者写成更加简洁的形式，提出 \[LeftAngleBracket]R\[VerticalSeparator]Subscript[\[Psi], nj]\[RightAngleBracket]，得到

\[\left(\nabla _n{}^2\text{ + }E_j\text{(R) + $\langle $}\psi _j\text{(R)$|$r$\rangle $}\nabla _n{}^2\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $ + $\langle $}\psi _j\text{(R)$|$r$\rangle \langle $r,R$|$}\hat{H}_{\text{mp}}\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $)$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ = }E_{\text{tot}}\text{$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ (Adiabatic Approximation)}\right.\]

如果我们进一步忽略质量极化，则有

\[\left(\nabla _n{}^2\text{ + }E_j\text{(R) + $\langle $}\psi _j\text{(R)$|$r$\rangle $}\nabla _n{}^2\text{$\langle $r$|$}\psi _j\text{(R)$\rangle $)$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ = }E_{\text{tot}}\text{$\langle $R$|$}\psi _{\text{nj}}\right\rangle\]

\[LeftAngleBracket]Subscript[\[Psi], j](R)\[VerticalSeparator]r\[RightAngleBracket]Subscript[\[Del], n]^2\[LeftAngleBracket]r\[VerticalSeparator]Subscript[\[Psi], j](R)\[RightAngleBracket] 被称为对角修正 (diagonal correction)，其数量级与 Subscript[E, j](R) 相比大致与电子-核质量比相同。

我们进一步将对角修正也忽略，则得到所谓的波恩-奥本海默近似(Born\[Dash]Oppenheimer approximation)

\[\left(\nabla _n{}^2\text{ + }E_j\text{(R))$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ = }E_{\text{tot}}\text{$\langle $R$|$}\psi _{\text{nj}}\text{$\rangle $ (Born$–$Oppenheimer approximation)}\right.\]

注意到 Subscript[\[Del], n]^2 = \[LeftAngleBracket]r,R\[VerticalSeparator]Subscript[Overscript[T, ^], n] = \[LeftAngleBracket]R\[VerticalSeparator]Subscript[Overscript[T, ^], n]，从而

\[\left(\hat{T}_n\text{ + }E_j\text{(R))$|$}\psi _{\text{nj}}\text{$\rangle $ = }E_{\text{tot}}\right|\psi _{\text{nj}}\rangle\]

显然，我们从形式上可以得到，此薛定谔方程描述的是原子核在一个 V = Subscript[E, j](R) 的势能中运动。该势能是电子薛定谔方程的解，被称为电子势能面(potential energy surface, PES)。解此方程，则得到近似的分子振动、转动能级等，进而可用于计算分子的光谱等。

显然，我们可以马上看出在BO近似下，PES与核质量无关，而在绝热近似下，PES由于对角修正的影响，将部分地与核质量有关。

所谓的 diagonal Born\[Dash]Oppenheimer correction (DBOC)，就是重新将之前忽略的 \[LeftAngleBracket]Subscript[\[Psi], j](R)\[VerticalSeparator]Subscript[\[Del], n]^2\[VerticalSeparator]Subscript[\[Psi], j](R)\[RightAngleBracket] 项引入。

BO近似在大多数情况下都是非常好的近似，除非电子薛定谔方程出现能量上非常相近的两个解。在这种情况下，将会有非绝热的效应出现。非绝热是目前一个非常有趣的问题，有大量的化学现象需要由非绝热过程完成。不过，考虑非绝热的情况会比较复杂，我们之后再来讨论相关的问题。现在，我们先从BO或者绝热近似开始，继续我们的讨论。

在继续之前，我们必须记住非常重要的一点：当没有BO近似的时候，所谓的PES是没有定义的，分子的几何形状，例如键长等，也都是没有精确定义的。

在BO近似的基础上，我们需要讨论究竟如何得到电子波函数。对于分子这样的多电子体系，薛定谔方程难以解析求解。由于忽略相对论性效应，我们只能先随意地引入电子自旋的概念：我们直接接受电子的自旋是1/2，有两种自旋状态，\[Alpha]或者\[Beta]，并且满足正交或者归一关系：

\[\text{$\langle \alpha |\beta \rangle $=}\delta _{\alpha \beta }\]

等到我们引入考虑相对论效应的狄拉克方程，自旋的概念就会更自然地浮现，不过现在我们就姑且把它当作一个既成的事实来使用，这不影响我们接下来的计算。

我们重新声明一下现在要进行的计算：在BO近似下，我们想要得到的是在核坐标固定的体系下的多电子运动方程，即其波函数和能级。因此，以下当我们使用BO近似的时候，我们不再写出R，比如我们直接记 Subscript[Overscript[H, ^], e](R) 为 Subscript[Overscript[H, ^], e]

\[E_i\text{=$\langle $}\psi _i\left|\hat{H}_e\right|\psi _i\rangle\]

解决这个问题由许多种方法，我们首先考虑变分法。

\[\text{== Variational Method ==}\]

定态薛定谔方程为

\[\hat{H}\text{$|\psi \rangle $=E$|\psi \rangle $}\]

当给定 Overscript[H, ^] 时，存在可数或不可数多个解，使得

\[\hat{H}\left|\psi _{\alpha }\text{$\rangle $=}E_{\alpha }\right|\psi _{\alpha }\text{$\rangle $, $\alpha $ = 0, 1, 2…}\]

且

\[E_0\text{ $<$ }E_1\text{ $<$ }E_2\text{ $<$ …}\]

为了简单起见，我们这里考虑可数多个解的情况，即能量是离散谱。对于连续谱，这里的证明亦可推广。

对于任意波函数\[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket]，Overscript[H, ^] 关于其的期望值中插入正交基，

\[\left\langle \tilde{\psi }\right|\hat{H}|\tilde{\psi }\text{$\rangle $= $\langle $}\tilde{\psi }|\hat{H}\sum _{\alpha =1}^{\infty } |\psi _{\alpha }\rangle \langle \psi _{\alpha }|\tilde{\psi }\text{$\rangle $=}\sum _{\alpha =1}^{\infty } \langle \tilde{\psi }|\hat{H}|\psi _{\alpha }\rangle \langle \psi _{\alpha }|\tilde{\psi }\text{$\rangle $=}\sum _{\alpha =1}^{\infty } E_{\alpha }\langle \tilde{\psi }|\psi _{\alpha }\rangle \langle \psi _{\alpha }|\tilde{\psi }\text{$\rangle $=}\sum _{\alpha =1}^{\infty } E_{\alpha }|\langle \psi _{\alpha }|\tilde{\psi }\rangle |^2\text{ $\geq $ }E_0\sum _{\alpha =1}^{\infty } |\langle \psi _{\alpha }|\tilde{\psi }\rangle |^2\text{ $\geq $ }E_0\]

等号成立当且仅当 \[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket] = \[VerticalSeparator]Subscript[\[Psi], 0]\[RightAngleBracket]。此即所谓变分原理。

注意到一般来说，试探矢量\[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket]的维数是无穷的。因此，直接优化\[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket]不是非常容易。为了克服这个问题，尤其是为了与现代计算机的结构相适应，我们一般采用线性变分法，即采用任意的一组正交归一完备基 {\[VerticalSeparator]i\[RightAngleBracket]}，它可以直接定义，也可以采用任何一个 Hamiltonian 的本征矢集。从而有

\[\left|\tilde{\psi }\text{$\rangle $ = }\sum _{i=1}^{\infty } \text{$|$i$\rangle \langle $i$|$}\tilde{\psi }\text{$\rangle $ = }\sum _{i=1}^{\infty } c_i\text{$|$i$\rangle $}\right.\]

而

\[\left\langle \tilde{\psi }\right|\hat{H}|\tilde{\psi }\text{$\rangle $ = $\langle $}\tilde{\psi }|\sum _{i=1}^{\infty } \text{$|$i$\rangle \langle $i$|$}\hat{H}\sum _{j=1}^{\infty } \text{$|$j$\rangle \langle $j$|$}\tilde{\psi }\text{$\rangle $=}\sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \langle \tilde{\psi }\text{$|$i$\rangle \langle $i$|$}\hat{H}\text{$|$j$\rangle \langle $j$|$}\tilde{\psi }\rangle\]

注意，\[VerticalSeparator]j\[RightAngleBracket] 并不是 Overscript[H, ^] 的本征态。上述问题现在仍是无穷的。不过，如果我们的基组选择非常不错(这就是我们为什么常常用其它体系的本征矢作为基组)，并且我们体系的能量不是太高，那么我们也许可以期望 \[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket] 展开在前几十或者几百个基矢上，就已经能概括\[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket]的大多数性质，从而在之后进行截断不会带来太大的误差。

即

\[\left\langle \tilde{\psi }\right|\hat{H}|\tilde{\psi }\text{$\rangle $ $\approx $ $\langle $}\tilde{\psi }\text{‘$|$}\hat{H}|\tilde{\psi }\text{‘$\rangle $ = $\langle $}\tilde{\psi }|\sum _{i=1}^N \text{$|$i$\rangle \langle $i$|$}\hat{H}\sum _{j=1}^N \text{$|$j$\rangle \langle $j$|$}\tilde{\psi }\text{$\rangle $ = }\sum _{i=1}^N \sum _{j=1}^N \langle \tilde{\psi }\text{$|$i$\rangle \langle $i$|$}\hat{H}\text{$|$j$\rangle \langle $j$|$}\tilde{\psi }\rangle\]

显然，当 N->\[Infinity]，\[LeftAngleBracket]Overscript[\[Psi], ~]’\[VerticalSeparator]Overscript[H, ^]\[VerticalSeparator]Overscript[\[Psi], ~]’\[RightAngleBracket]->\[LeftAngleBracket]Overscript[\[Psi], ~]\[VerticalSeparator]Overscript[H, ^]\[VerticalSeparator]Overscript[\[Psi], ~]\[RightAngleBracket]。

定义矩阵 H，Subscript[H, ij]=\[LeftAngleBracket]i\[VerticalSeparator]Overscript[H, ^]\[VerticalSeparator]j\[RightAngleBracket]，矢量Overscript[c, \[RightVector]]，Subscript[Overscript[c, \[RightVector]], i]=Subscript[c, i]，则

\[\left\langle \tilde{\psi }\text{‘$|$}\hat{H}\right|\tilde{\psi }\text{‘$\rangle $ = }\sum _{i=1}^N \sum _{j=1}^N \langle \tilde{\psi }\text{$|$i$\rangle \langle $i$|$}\hat{H}\text{$|$j$\rangle \langle $j$|$}\tilde{\psi }\text{$\rangle $ = }\sum _{i=1}^N \sum _{j=1}^N c_iH_{\text{ij}}c_j\text{= }\overset{\rightharpoonup }{c}^{\dagger }H\overset{\rightharpoonup }{c}\]

对\[VerticalSeparator]Overscript[\[Psi], ~]’\[RightAngleBracket]变分，由于存在约束

\[\left\langle \tilde{\psi }\text{‘$|$}\tilde{\psi }\text{‘$\rangle $=1}\right.\]

因此采用 Lagrange’s method of undetermined multipliers，

\[\text{L(}c_1\text{, }c_2\text{, }c_3\text{, …, }c_N\text{, E) = $\langle $}\tilde{\psi }\text{‘$|$}\hat{H}|\tilde{\psi }\text{‘$\rangle $ – E$\cdot $($\langle $}\tilde{\psi }\text{‘$|$}\tilde{\psi }\text{‘$\rangle $ – 1) = }\sum _{i=1}^N \sum _{j=1}^N c_iH_{\text{ij}}c_j\text{ – E(}\sum _{i=1}^N c_i{}^2\text{-1)}\]

利用 Subscript[H, ij]=Subscript[H, ji]，

\[\frac{\partial L}{\partial c_k}=\frac{\partial }{\partial c_k}\left(\sum _{i=1}^N \sum _{j=1}^N c_iH_{\text{ij}}c_j\text{ – E(}\sum _{i=1}^N c_i{}^2\text{-1))=}\frac{\partial }{\partial c_k}\right(\sum _{i=1}^N \sum _{j=1}^N c_iH_{\text{ij}}c_j\text{)-E}\frac{\partial }{\partial c_k}(\sum _{i=1}^N c_i{}^2\text{-1)=}\sum _{i=1}^N c_iH_{\text{ik}}\text{ + }\sum _{j=1}^N H_{\text{kj}}c_j\text{ – 2}\text{Ec}_k\text{ = }\sum _{i=1}^N \text{ 2}c_iH_{\text{ik}}\text{ – 2}\text{Ec}_k\text{ = 0}\]

\[\sum _{i=1}^N H_{\text{ik}}c_k\text{ – }\text{Ec}_k\text{ = 0}\]

显然

\[\sum _{k=1}^N \sum _{i=1}^N H_{\text{ik}}c_k\text{ – }\sum _{k=1}^N \text{Ec}_k\text{ = 0}\]

即

\[H\overset{\rightharpoonup }{c}\text{ = E}\overset{\rightharpoonup }{c}\]

此为矩阵本征值问题，即我们需要对角化矩阵 H

\[U^{\dagger }\text{HU = E}\]

记

\[\text{C = U}\]

则

\[\text{HC = CE}\]

此即变分法所得的结果。对应的最小本征值即为变分法的最低能量，该本征值对应的本征矢就是基态波函数展开在基组上的组合系数。

如果我们并未采用正交的基组，则我们还必须在结果中加入重叠积分项 Subscript[S, ij]，即此时方程变为

\[H\overset{\rightharpoonup }{c}\text{ = ES}\overset{\rightharpoonup }{c}\]

所得结果为

\[\text{HC = SCE}\]

从以上的分析我们可以看出，变分法求波函数的主要计算花费在计算 \[LeftAngleBracket]i\[VerticalSeparator]Overscript[H, ^]\[VerticalSeparator]j\[RightAngleBracket] 这个积分上。一旦矩阵 H 构建完成，对角化用现代的计算机非常快速。

\[\text{[CODE EXAMPLE 1]}\]

\[\text{[In this example, we demonstrate the linear variational technique that solves arbitrary 1-dimensional systems numerically using eigenfunctions of particle-in-a-box or eigenfunctions of harmonic oscillator as a basis set]}\]

\[\text{===Self-Consistent Field===}\]

变分原理告诉我们，当我们使用一个满足对称性要求的试探波函数代替\[VerticalSeparator]Subscript[Overscript[\[Psi], ~], i]\[RightAngleBracket]来计算

\[\tilde{E}\text{=$\langle $}\tilde{\psi }_i\left|\hat{H}_e\right|\tilde{\psi }_i\rangle\]

时，使 Overscript[E, ~] 最小的试探波函数 \[VerticalSeparator]Subscript[Overscript[\[Psi], ~], i]\[RightAngleBracket] 就是真实的基态波函数 \[VerticalSeparator]Subscript[\[Psi], i]\[RightAngleBracket]。因此，我们来考虑怎样构造一个N电子波函数。显然，N电子波函数是一个4N自由度的函数。我们不妨这样来构造：

假如体系中的电子均独立运动，我们有单电子近似哈密顿算符

\[\hat{H}_e\text{ = }\sum _{i=1}^N \hat{h}_i\]

从而若解

\[\hat{h}_i\left|\chi _i\text{$\rangle $ = }\epsilon _i\right|\chi _i\rangle\]

可以得到无穷多个单电子波函数 Subscript[\[Chi], i](Subscript[x, j])，称为自旋轨道，若将N个自旋轨道相乘，则得到一个4N自由度的函数，称为 Hartree Product

\[\Psi _{\text{HP}}=\chi _i\left(x_1\right)\chi _j(x_2\text{)…}\chi _k(x_N)\]

所谓的对称性要求，实际上就是变分域的约束条件，即电子作为费米子，满足波函数对交换电子反对称。

为了满足对称性要求，我们的波函数需要采取一定的特殊形式。一种从任意形式构造反对称性的方法是使用 Slater determinant。举个例子：如果我们用N个正交的自旋轨道和N个电子，那么我们可以这样构造一个反对称的电子波函数：

\[\Phi _{\text{SD}}=\frac{1}{\sqrt{N!}}\text{$|\backslash $!$\backslash $($\backslash $*GridBox[$\{\{$RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}i{”}], RowBox[$\{$(, 1, )$\}$]$\}$], RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}j{”}], RowBox[$\{$(, 1, )$\}$]$\}$], …, RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}k{”}], RowBox[$\{$(, 1, )$\}$]$\}$]$\}$,$\{$RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}i{”}], RowBox[$\{$(, 2, )$\}$]$\}$], RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}j{”}], RowBox[$\{$(, 2, )$\}$]$\}$], …, RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}k{”}], RowBox[$\{$(, 2, )$\}$]$\}$]$\}$,$\{$…, …, …, …$\}$,$\{$RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}i{”}], RowBox[$\{$(, N, )$\}$]$\}$], RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}j{”}], RowBox[$\{$(, N, )$\}$]$\}$], …, RowBox[$\{$SubscriptBox[{“}$\chi ${”}, {“}k{”}], RowBox[$\{$(, N, )$\}$]$\}$]$\}$}\]

\[\text{ $\}$]$\backslash $)$|$}\]

简记为

\[\left|\chi _1\chi _2\text{…}\chi _N\right\rangle\]

Slater determinant 具有如下的性质：

\[\text{$|$…}\chi _n\text{…}\chi _m\text{…$\rangle $ = -$|$…}\chi _m\text{…}\chi _n\text{…$\rangle $}\]

\[\text{$\langle $…}\chi _m\text{…$|$…}\chi _n\text{…$\rangle $ = }\delta _{\text{nm}}\]

显然，任意多个 Slater determinant 的线性组合也都是满足反对称的。

由于真实波函数\[VerticalSeparator]Subscript[\[Psi], i]\[RightAngleBracket]的自由度为4N，既然自旋轨道是自由度为4的完备基组，那么真实波函数\[VerticalSeparator]Subscript[\[Psi], i]\[RightAngleBracket]应该可以展开在N个自旋轨道的乘积上

\[\psi _i\left(x_1\text{, }x_2\text{, …, }x_N\text{)=}\sum _{i,j,\text{…},k=1}^{\infty } b_{\text{ij}\text{…}k}\chi _i\right(x_1)\chi _j(x_2\text{)…}\chi _k(x_N)\]

Slater determinant是4N自由度的函数，显然，Slater determinant对于4N自由度任意波函数并不完备，不过，我们可以证明Slater determinant对于4N自由度的反对称波函数是完备的。因此，真实波函数\[VerticalSeparator]Subscript[\[Psi], i]\[RightAngleBracket]同样可以展开在Slater determinant上。

如果我们允许变分域为全反对称域的话，那么我们通过变分法就能得到准确的电子波函数，这个办法被叫做 Full CI。不过遗憾的是，以后我们会看到，这个想法可能

会导致我们的计算非常困难。

如果我们进一步限制变分域，那么显然通过变分法，只能得到能量的一个上界。如果用这个上界来估计体系的真实电子波函数能量，那么实际上我们就又做了另一个近似。不过，现在先让我们做一个非常大胆的限制，即把变分域限制为仅仅由一个上述的 slater determinant 构成的波函数，从而我们只能改变N个上述的自旋轨道。这实际上意味着，电子之间的所有耦合作用都被忽略了，即电子之间的排斥作用仅仅被视作一种平均的作用。这一近似相当于展开在Slater determinant上并且只取最重要的一项，虽然可以大大简化计算，但它是非常不精确的近似，因此我们以后还要转回来看怎样改进它。

那么，我们已经做了如此多的近似，是否可以开始计算了呢？非常遗憾，即使对于仅仅一个Slater determinant，我们也还是无法处理，因为组成Slater determinant的自旋轨道有无穷多的选择方式。如前面变分法的介绍所述，为了变分，我们还需要对基组进行截断。通常来说，我们取K个基组，从中取N个组成我们的Slater determinant进行变分，将得到K个正交的自旋轨道Subscript[\[Chi], 1], Subscript[\[Chi], 2], …Subscript[\[Chi], K]。,而\[VerticalSeparator]Subscript[\[Chi], 1]Subscript[\[Chi], 2]…Subscript[\[Chi], N]\[RightAngleBracket]就是能量最低的Slater determinant，我们将其作为近似电子基态函数。当 K->\[Infinity] 时，得到的近似结果，被称为 Hartree-Fock Limit。

从以上的分析我们可以看出，在 Hartree-Fock Limit 和 Full CI 的条件下，我们就能通过变分法得到准确的波函数。不过为了能进行实际的计算，我们进行了大量近似。

在这些近似下，我们可以现在来推导 Hartree-Fock 方程了。Hartree-Fock 方程是通过K个基组计算K个自旋轨道的一种方法，由于Fock算符中具有对体系的自指，该方程组一般用迭代法求解，即所谓的Self-Consistent Field (SCF)。不过，在推导 Hartree-Fock 方程之前，为了方便，我们先看看如何简化 Slater determinant 的能量期望计算

\[E_{\text{SD}}\text{=$\langle $}\Phi _{\text{SD}}\left|\hat{H}_e\right|\Phi _{\text{SD}}\rangle\]

\[\hat{H}_e\text{ = }\hat{T}_e\text{ + }\hat{V}_{\text{ee}}\text{ + }\hat{V}_{\text{ne}}\text{ + }\hat{V}_{\text{nn}}\]

\[\text{$\langle $r,R$|$}\hat{T}_e\text{ = -}\sum _{i=1}^{N_{\text{elec}}} \frac{1}{2}\nabla _i{}^2\]

\[\text{$\langle $r,R$|$}\hat{V}_{\text{ee}}\text{ = }\sum _{i=1}^{N_{\text{elec}}} \sum _{j>i}^{N_{\text{elec}}} \frac{1}{\left|r_i-r_j\right|}\]

\[\text{$\langle $r,R$|$}\hat{V}_{\text{ne}}\text{ = -}\sum _{i=1}^{N_{\text{elec}}} \sum _{a=1}^{N_{\text{nuc}}} \frac{Z_a}{\left|R_a-r_i\right|}\]

\[\text{$\langle $r,R$|$}\hat{V}_{\text{nn}}\text{ = }\sum _{a=1}^{N_{\text{nuc}}} \sum _{b>a}^{N_{\text{nuc}}} \frac{Z_aZ_b}{\left|R_a-R_b\right|}\]

由于R在此时是参量，以上算符可以重新整理为单电子和双电子算符

\[\text{$\langle $r,R$|$}\hat{h}_i\text{ = -}\frac{1}{2}\nabla _i{}^2\text{ – }\sum _{a=1}^{N_{\text{nuc}}} \frac{Z_a}{\left|R_a-r_i\right|}\]

\[\text{$\langle $r,R$|$}\hat{g}_{\text{ij}}\text{ = }\frac{1}{\left|r_i-r_j\right|}\]

从而

\[\hat{H}_e\text{ = }\sum _{i=1}^{N_{\text{elec}}} \hat{h}_i\text{ + }\sum _{j>i}^{N_{\text{elec}}} \hat{g}_{\text{ij}}\text{ + }\hat{V}_{\text{nn}}\text{ = }\hat{h}\text{ + }\hat{g}\text{ + }\hat{V}_{\text{nn}}\]

如之前变分法所述，为了变分，我们需要构造矩阵 Subscript[H, ijk…lpqr…s]，我们引入简记

\[\left\langle \Phi _{\text{SD1}}\right|\hat{H}_e|\Phi _{\text{SD2}}\text{$\rangle $ = $\langle $}\chi _i\chi _j\text{…}\chi _k|\hat{H}_e|\chi _p\chi _q\text{…}\chi _r\text{$\rangle $ = $\langle $ij…k$|$}\hat{H}_e\text{$|$pq…r$\rangle $}\]

注意到

\[\text{$\langle $mn…o$|$}\hat{H}_e\text{$|$pq…r$\rangle $ = $\langle $mn…o$|$}\sum _{i=1}^{N_{\text{elec}}} \hat{h}_i\text{ + }\sum _{j>i}^{N_{\text{elec}}} \hat{g}_{\text{ij}}\text{ + }\hat{V}_{\text{nn}}\text{$|$pq…r$\rangle $ = }\sum _{i=1}^{N_{\text{elec}}} \text{$\langle $mn…o$|$}\hat{h}_i\text{$|$pq…r$\rangle $ + }\sum _{j>i}^{N_{\text{elec}}} \text{$\langle $mn…o$|$}\hat{g}_{\text{ij}}\text{$|$pq…r$\rangle $ + $\langle $mn…o$|$}\hat{V}_{\text{nn}}\text{$|$pq…r$\rangle $}\]

即我们需要计算

\[\text{$\langle $mn…o$|$}\hat{h}_i\text{$|$pq…r$\rangle $}\]

\[\text{$\langle $mn…o$|$}\hat{g}_{\text{ij}}\text{$|$pq…r$\rangle $}\]

\[\text{$\langle $mn…o$|$}\hat{V}_{\text{nn}}\text{$|$pq…r$\rangle $}\]

首先引入一些简记：

对自旋轨道，

\[\text{$\langle $m$|$}\hat{h}_i\text{$|$n$\rangle $ = [m$|$}\hat{h}_i\text{$|$n] = $\int $}\chi _m\left(x_i\right)\hat{h}_i(r_i)\chi _n(x_i\text{)d}x_i\]

\[\text{$\langle $mn$|$pq$\rangle $ = [mp$|$nq] = $\int \int $}\chi _m\left(x_i\right)\chi _n(x_j)\hat{g}_{\text{ij}}(r_i\text{, }r_j)\chi _p(x_i)\chi _q(x_j\text{)d}x_idx_j\text{ (Note that the integration is independent of i, j (they are pseudo indices), so we can also just use 1, 2 instead of i, j)}\]

\[\text{$\langle $mn$||$pq$\rangle $ = $\langle $mn$|$pq$\rangle $ – $\langle $mn$|$qp$\rangle $}\]

对空间轨道，

\[\text{(m$|$}\hat{h}_i\text{$|$n) = $\int $}\psi _m\left(r_i\right)\hat{h}_i(r_i)\psi _n(r_i\text{)d}r_i\]

\[\text{(mn$|$pq) = $\int \int $}\psi _m\left(r_i\right)\psi _n(r_i)\hat{g}_{\text{ij}}(r_i\text{, }r_j)\psi _p(r_j)\psi _q(r_j\text{)d}r_idr_j\]

\[\text{(mn$||$pq) = (mn$|$pq) – (mn$|$qp)}\]

\[J_{\text{ij}}\text{=(ii$|$jj) (Coulumb integrals)}\]

\[K_{\text{ij}}\text{=(ij$|$ji) (Exchange integrals)}\]

我们现在来计算能量

对于单电子算符，

\[\text{$\langle $…mn…$|$}\hat{h}_i\text{$|$…mn…$\rangle $ = $\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $}\]

\[\text{$\langle $…mn…$|$}\hat{h}_i\text{$|$…pn…$\rangle $ = 0}\]

从而

\[\text{$\langle $…mn…$|$}\hat{h}\text{$|$…mn…$\rangle $ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $}\]

\[\text{$\langle $…mn…$|$}\hat{h}\text{$|$…pn…$\rangle $ = $\langle $m$|$}\hat{h}\text{$|$p$\rangle $}\]

\[\text{$\langle $…mn…$|$}\hat{h}\text{$|$…pq…$\rangle $ = 0}\]

对于双电子算符，

\[\text{$\langle $…ij…$|$}\hat{g}_{\text{ij}}\text{$|$…ij…$\rangle $ = $\langle $ij$|$ij$\rangle $}\]

\[\text{$\langle $…ij…$|$}\hat{g}_{\text{ij}}\text{$|$…ji…$\rangle $ = $\langle $ij$|$ji$\rangle $ }\]

\[\text{$\langle $…k…$|$}\hat{g}_{\text{ij}}\text{$|$…p…$\rangle $ = 0}\]

从而

\[\text{$\langle $…mn…$|$}\hat{g}\text{$|$…mn…$\rangle $ =}\frac{1}{2}\sum _{p=1}^N \sum _{q=1}^N \text{$\langle $pq$||$pq$\rangle $}\]

\[\text{$\langle $…mn…$|$}\hat{g}\text{$|$…pn…$\rangle $ = }\sum _{q=1}^N \text{$\langle $mq$||$pq$\rangle $}\]

\[\text{$\langle $…mn…$|$}\hat{g}\text{$|$…pq…$\rangle $ = $\langle $mn$||$pq$\rangle $}\]

注意，以上结果要求先将determinant的序号排列为严格升序。比如说 \[VerticalSeparator]1423\[RightAngleBracket] 应该先经过重新排序

\[\text{$|$1423$\rangle $ = -$|$1324$\rangle $ = $|$1234$\rangle $}\]

再套用以上规则。

一般，我们用 a, b, … 代表占据轨道， r, s, … 代表虚轨道。非基态的 Slater determinant 也可以看作是将电子从占据轨道上 “激发” 到虚轨道得来，从而可记作

\[\left|\Phi _0\text{$\rangle $ = $|$12…ab…N$\rangle $}\right.\]

\[\left|\Phi _{\text{ab}}^{\text{rs}}\text{$\rangle $ = $|$12…rs…N$\rangle $}\right.\]

用之前的结果我们可以比较容易地计算这些 Slater determinant 的各种积分。

现在，我们先考虑\[VerticalSeparator]Subscript[\[CapitalPhi], 0]\[RightAngleBracket] = \[VerticalSeparator]12…ab…N\[RightAngleBracket] 的 Hartree-Fock 能量。

\[E_0\text{ = $\langle $12…ab…N$|$}\hat{H}_e\text{$|$12…ab…N$\rangle $ = $\langle $12…ab…N$|$}\hat{h}\text{ + }\hat{g}\text{ + }\hat{V}_{\text{nn}}\text{$|$12…ab…N$\rangle $ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $ + }\sum _{i=1}^N \sum _{j>i}^N \text{($\langle $ij$|$ij$\rangle $ – $\langle $ij$|$ji$\rangle $) + }V_{\text{nn}}\text{ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $ +}\frac{1}{2}\sum _{i,j=1}^N \text{($\langle $ij$|$ij$\rangle $ – $\langle $ij$|$ji$\rangle $) + }V_{\text{nn}}\]

如果我们进一步限制体系是一个 restricted closed shell 体系，即2K个自旋轨道是K个空间轨道直接乘以不同的自旋函数得到，并且占据轨道都按自旋成对，那么可以进一步简化 Subscript[E, 0] 为

\[E_0\text{ = 2}\sum _{i=1}^{N/2} \text{(i$|$}\hat{h}_i\text{$|$i) + }\sum _{i,j=1}^{N/2} \text{(2(ii$|$jj) – (ij$|$ji)) + }V_{\text{nn}}\text{ = 2}\sum _{a=1}^{N/2} h_{\text{aa}}\text{ + }\sum _{a,b=1}^{N/2} 2J_{\text{ab}}\text{ – }K_{\text{ab}}\]

我们可以这样来帮助记忆以上的式子：空间轨道 Subscript[\[Psi], i] 上的每个电子贡献一个动能项 Subscript[h, ii]，每对电子，假如它们分别处在 Subscript[\[Psi], i]，Subscript[\[Psi], j] 上，贡献一个库伦积分项 Subscript[J, ij]，每对自旋相同的电子，贡献一个交换积分项 Subscript[K, ij]。由于自旋相同的电子是总电子数的一半，故交换积分项比库伦积分项少一半。

必须注意：以上的说法只是为了帮助记忆，所谓的”库伦”，”交换”，都只是人造的概念，其本质是运动耦合，电子间并没有什么真正的”交换”作用。

任意的非基态的 Slater determinant 的能量，也可以用类似以上的流程计算。

(TODO: 简单介绍二次量子化，以及用二次量子化表述再推导一遍以上的结果)

== 自旋问题 ==

我们注意到，对于任意的 Slater determinant，它可能是总角动量平方算符 Overscript[J^2, ^] 的本征态，也有可能不是。对于 Restricted Slater Determinant，我们可以通过线性组合非本征态的方法得到本征态，即所谓的 spin-adapted configuration。

例如，对于 Minimal Basis Subscript[H, 2] ，\[VerticalSeparator]Subscript[\[CapitalPhi], 0]\[RightAngleBracket]，\[VerticalSeparator]Subsuperscript[\[CapitalPhi], 1Overscript[1, _], 2Overscript[2, _]]\[RightAngleBracket]，\[VerticalSeparator]Subsuperscript[\[CapitalPhi], 1, Overscript[2, _]]\[RightAngleBracket]，\[VerticalSeparator]Subsuperscript[\[CapitalPhi], Overscript[1, _], 2]\[RightAngleBracket] 是Overscript[J^2, ^] 的本征态，而\[VerticalSeparator]Subsuperscript[\[CapitalPhi], 1, 2]\[RightAngleBracket]，\[VerticalSeparator]Subsuperscript[\[CapitalPhi], Overscript[1, _], Overscript[2, _]]\[RightAngleBracket]则不是。

\[\left|\Phi _1^2\text{$\rangle $ = $\backslash $!$\backslash $($\backslash $*GridBox[$\{\{$UnderscriptBox[{“}$\uparrow ${”}, {“}$\_${”}]$\}$,$\{$UnderscriptBox[{“}$\downarrow ${”}, {“}$\_${”}]$\}$}\right.\]

\[\text{ $\}$]$\backslash $)}\]

\[\left|\Phi _{\bar{1}}^{\bar{2}}\text{$\rangle $ = $\backslash $!$\backslash $($\backslash $*GridBox[$\{\{$UnderscriptBox[{“}$\downarrow ${”}, {“}$\_${”}]$\}$,$\{$UnderscriptBox[{“}$\uparrow ${”}, {“}$\_${”}]$\}$}\right.\]

\[\text{ $\}$]$\backslash $)}\]

不过，通过其线性组合，可以得到 Overscript[J^2, ^] 的本征态(与反对称化做法相同)

Unrestricted Slater Determinant 一般不是 Overscript[J^2, ^] 的本征态，并且无法通过线性组合一小部分 determinant 来形成 spin-adapted configuration。相比自旋纯态，UHF基态总是被更高多重度的自旋态污染，因此，Overscript[J^2, ^] 的期望值总是偏大的。

\[\text{=== Hartree-Fock ===}\]

在许多初级的化学课本上，分子中的电子各自占据不同的分子轨道，这实际上是一个近似，叫做Hartree-Fock近似，或者叫分子轨道近似。现在，我们可以暂时就把Hartree-Fock theory 看作是 Single determinant theory。

Koopmans’ theorem：单行列式近似下，自旋轨道能量可以理解为在该轨道失去或得到一个电子的近似的电离能或电子亲和能。

在得到或失去一个电子之后，体系的所有轨道波函数都应该发生变化，这叫做relaxation，因为该电子对其它电子的排斥作用改变了。而在Koopmans’ theorem中，波函数仿佛被冻结了。因此，用Koopmans’ theorem估算的电离能或亲和能会绝对值偏大。不过，由于单行列式近似总是抬高体系的总能量(证明见变分法章节)，因此对电离能的估算正负误差会互相抵消一部分，从而反而得到和实验值比较接近的结果，而对亲和能的误差则两者叠加，得到非常差的结果。

Brillouin’s theorem：单激发行列式不会直接与基态行列式作用，而是通过高次项耦合进来。因此，考虑对基态的修正时，单独引入单激发行列式没有作用，首先应该考虑的是双激发行列式(CID)，然后才考虑 CISD 等等

现在我们来变分。显然，任意给定的 Slater determinant 由选定的N个自旋轨道波函数完全确定。因此，变分 Slater determinant 等价于寻找N个自旋轨道使基态determinant能量最低。在之前，我们已经得到 determinant 能量的表达式

\[E_0\text{ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $ +}\frac{1}{2}\sum _{i,j=1}^N \text{($\langle $ij$|$ij$\rangle $ – $\langle $ij$|$ji$\rangle $) + }V_{\text{nn}}\]

我们可以将其表达为一个更方便的形式，引入库伦算符 (Coulomb Operator) Subscript[Overscript[J, ^], i] 和交换算符 (Exchange Operator) Subscript[Overscript[K, ^], i]，

\[\hat{J}_i\text{$|$j(}x_m\text{)$\rangle $ = $\langle $i(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$i(}x_n\text{)$\rangle |$j(}x_m\text{)$\rangle $ (Note: like the expression in $\langle $ij$|$ij$\rangle $, m $\&$ n are just psuedo indices)}\]

\[\hat{K}_i\text{$|$j(}x_m\text{)$\rangle $ = $\langle $i(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$j(}x_n\text{)$\rangle |$i(}x_m\text{)$\rangle $}\]

用以上的定义，有

\[\text{$\langle $j$|$}\hat{J}_i\text{$|$j$\rangle $ = $\langle $j(}x_m\text{)$|\langle $i(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$i(}x_n\text{)$\rangle |$j(}x_m\text{)$\rangle $ = $\langle $ji$|$ji$\rangle $}\]

\[\text{$\langle $i$|$}\hat{J}_j\text{$|$i$\rangle $ = $\langle $i(}x_m\text{)$|\langle $j(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$j(}x_n\text{)$\rangle |$i(}x_m\text{)$\rangle $ = $\langle $ij$|$ij$\rangle $}\]

\[\langle \alpha |\hat{J}_j\text{$|\beta \rangle $ = $\langle \alpha $(}x_m\text{)$|\langle $j(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$j(}x_n\text{)$\rangle |\beta $(}x_m\text{)$\rangle $ = $\langle \alpha $j$|\beta $j$\rangle $}\]

\[\text{$\langle $j$|$}\hat{K}_i\text{$|$j$\rangle $ = $\langle $j(}x_m\text{)$|\langle $i(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$j(}x_n\text{)$\rangle |$i(}x_m\text{)$\rangle $ = $\langle $ji$|$ij$\rangle $}\]

\[\text{$\langle $i$|$}\hat{K}_j\text{$|$i$\rangle $ = $\langle $i(}x_m\text{)$|\langle $j(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|$i(}x_n\text{)$\rangle |$j(}x_m\text{)$\rangle $ = $\langle $ij$|$ji$\rangle $}\]

\[\langle \alpha |\hat{K}_j\text{$|\beta \rangle $ = $\langle \alpha $(}x_m\text{)$|\langle $j(}x_n\text{)$|$}\hat{g}_{\text{mn}}\text{$|\beta $(}x_n\text{)$\rangle |$j(}x_m\text{)$\rangle $ = $\langle \alpha $j$|$j$\beta \rangle $}\]

(这里选了ij作为标记真是太失败了。。啥都看不清[Facepalm]，但我不想改了)

从而能量可以写为

\[E_0\text{ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $ +}\frac{1}{2}\sum _{i,j=1}^N \text{($\langle $i$|$}\hat{J}_j\text{$|$i$\rangle $ – $\langle $i$|$}\hat{K}_j\text{$|$i$\rangle $) + }V_{\text{nn}}\]

Subscript[V, nn] 在BO近似下显然是常数，因此它不影响我们的变分，我们最后再把它加回来就可以，现在先姑且忽略它

\[\tilde{E}_0\text{ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{h}_i\text{$|$i$\rangle $ +}\frac{1}{2}\sum _{i,j=1}^N \text{($\langle $i$|$}\hat{J}_j\text{$|$i$\rangle $ – $\langle $i$|$}\hat{K}_j\text{$|$i$\rangle $)}\]

定义 Fock 算符

\[\hat{F}_i\text{ = }\hat{h}_i\text{ + }\sum _{j=1}^N \text{ ( }\hat{J}_j\text{ – }\hat{K}_j\text{ )}\]

则能量可以写成

\[\tilde{E}_0\text{ = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{F}_i\text{$|$i$\rangle $}\]

现在我们终于可以来变分了，由于存在约束

\[\text{$\langle $i$|$j$\rangle $ = }\delta _{\text{ij}}\]

需要用拉格朗日乘因子法

\[\left.\text{L = }\tilde{E}_0\text{ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{($\langle $i$|$j$\rangle $ – }\delta _{\text{ij}}\text{) = }\sum _{i=1}^N \text{$\langle $i$|$}\hat{F}_i\text{$|$i$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{($\langle $i$|$j$\rangle $ – }\delta _{\text{ij}}\right)\]

\[\text{$\delta $L = }\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i$\rangle $ + }\sum _{i=1}^N \text{$\langle $i$|$}\hat{F}_i\text{$|\delta $i$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{$\langle \delta $i$|$j$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{$\langle $i$|\delta $j$\rangle $ = 0}\]

\[\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{$\langle \delta $i$|$j$\rangle $ + }\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i}\rangle ^*\text{ – }\sum _{i,j=1}^N \lambda _{\text{ji}}\text{$\langle $j$|\delta $i$\rangle $ = 0}\]

\[\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{$\langle \delta $i$|$j$\rangle $ + }\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i}\rangle ^*\text{ – }\sum _{i,j=1}^N \lambda _{\text{ji}}^*\text{$\langle \delta $i$|$j}\rangle ^*\text{ = 0}\]

变分 \[LeftAngleBracket]\[Delta]i\[VerticalSeparator] 或 \[LeftAngleBracket]\[Delta]i\[VerticalSeparator]^* 都应该使上式取0，故上式两部分必须各自取0

\[\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ij}}\text{$\langle \delta $i$|$j$\rangle $ = 0}\]

\[\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i}\rangle ^*\text{ – }\sum _{i,j=1}^N \lambda _{\text{ji}}\text{$\langle \delta $i$|$j}\rangle ^*\text{ = 0}\]

即

\[\hat{F}_i\text{$|$i$\rangle $ = }\sum _{j=1}^N \lambda _{\text{ij}}\text{$|$j$\rangle $}\]

此即 Hartree-Fock 方程。我们注意到 Hartree-Fock 方程具有类似本征值方程的形式，事实上，我们确实可以将其变成一个本征值方程的形式，叫做 Canonical Hartree-Fock 方程。我们现在来看怎么做这件事

首先，我们取刚才两式的第二式的共轭，得到

\[\sum _{i=1}^N \text{$\langle \delta $i$|$}\hat{F}_i\text{$|$i$\rangle $ – }\sum _{i,j=1}^N \lambda _{\text{ji}}^*\text{$\langle \delta $i$|$j$\rangle $ = 0}\]

将其与第一式相减，得到

\[\sum _{i,j=1}^N \left(\lambda _{\text{ij}}\text{ – }\lambda _{\text{ji}}^*\text{)$\langle \delta $i$|$j$\rangle $}\right.\]

\[\text{i.e.}\]

\[\lambda _{\text{ij}}\text{ = }\lambda _{\text{ji}}^*\]

即 Subscript[\[Lambda], ij] 是一个 rank-N Hermitian 矩阵 \[CapitalLambda] 的矩阵元。从而有使 \[CapitalLambda] 对角化的酉矩阵 U，s.t.

\[\text{E = }U^{\dagger }\text{$\Lambda $U}\]

\[\text{E = diag([}\varepsilon _1\text{, }\varepsilon _2\text{, }\varepsilon _3\text{, …, }\varepsilon _N\text{])}\]

\[\text{i.e.}\]

\[E_{\text{ij}}\text{ = }\varepsilon _i\delta _{\text{ij}}\]

对基组也进行对应的酉变换，从而 Hartree-Fock 方程化为

\[\hat{F}_i\text{$|$i$\rangle $ = }\sum _{j=1}^N E_{\text{ij}}\text{$|$j$\rangle $ = }\sum _{j=1}^N \varepsilon _i\delta _{\text{ij}}\text{$|$j$\rangle $ = }\varepsilon _i\text{$|$i$\rangle $}\]

简单地就写作

\[\hat{F}_i\text{$|$i$\rangle $ = }\varepsilon _i\text{$|$i$\rangle $}\]

此为 Canonical Hartree-Fock equation.

需要注意，Subscript[这里得到的\[CurlyEpsilon], i]是轨道能量，由于电子间存在耦合，

\[\tilde{E}_0\text{ $\neq $ }\sum _{i=1}^N \varepsilon _i\]

Subscript[Overscript[E, ~], 0] 需要用求得的自旋轨道构造的 Slater determinant 来计算。

从上述分析我们看到，满足 Canonical Hartree-Fock equation 的自旋轨道，就等价于变分 Slater determinant 求得的轨道。

既然如此，我们再来求 Canonical Hartree-Fock equation 的解。一种方法是再用一次变分法，得到所谓的 Roothaan 方程或者 Pople-Nesbet 方程。引入任意一个基组，如我们一开始介绍的变分法，本征值方程可化为

\[\text{FC = SCE}\]

其中，S 是基组的重叠矩阵。现在我们遇到了一个问题：前面求解方程 HC = SCE 时，Overscript[H, ^] 中并不含波函数的信息，因此我们可以直接 evaluate \[LeftAngleBracket]i\[VerticalSeparator]Overscript[H, ^]\[VerticalSeparator]j\[RightAngleBracket]

而 FC = SCE 虽然形式上与之相同，但 Overscript[F, ^] 的定义中包含对求解结果的自指，因此这是非常非线性的方程。一般来说，我们常常用迭代的办法求解，将上一次的计算结果作为猜测代入下一次计算，直到结果不再变化为止。这种方法被称为 self-consistent field (SCF) 方法

Fock 算符的定义是

\[\hat{F}_i\text{ = }\hat{h}_i\text{ + }\sum _{j=1}^N \text{ ( }\hat{J}_j\text{ – }\hat{K}_j\text{ )}\]

假设我们引入的 K 个基组为 {\[VerticalSeparator]\[Alpha]\[RightAngleBracket]}，则

\[S_{\alpha \beta }\text{ = $\langle \alpha |\beta \rangle $}\]

\[F_{\alpha \beta }\text{ = $\langle \alpha |$}\hat{F}_i\text{$|\beta \rangle $ = $\langle \alpha |$}\hat{h}_i\text{ + }\sum _{j=1}^N \text{ ( }\hat{J}_j\text{ – }\hat{K}_j\text{ )$|\beta \rangle $ = $\langle \alpha |$}\hat{h}_i\text{$|\beta \rangle $ + }\sum _{j=1}^N \text{ $\langle \alpha |$}\hat{J}_j\text{$|\beta \rangle $ – $\langle \alpha |$}\hat{K}_j\text{$|\beta \rangle $ = }H_{\alpha \beta }^{\text{core}}\text{ + }\sum _{j=1}^N \text{ $\langle \alpha $j$|\beta $j$\rangle $ – $\langle \alpha $j$|$j$\beta \rangle $ = }H_{\alpha \beta }^{\text{core}}\text{ + }\tilde{G}_{\alpha \beta }\]

注意到

\[H_{\alpha \beta }^{\text{core}}\text{ = $\langle \alpha |$}\hat{h}_i\text{$|\beta \rangle $ = }T_{\alpha \beta }\text{ + }V_{\alpha \beta }^{\text{ne}}\]

是与解无关的，故可以直接 evaluate，而计算

\[\tilde{G}_{\alpha \beta }\text{ = }\sum _{j=1}^N \text{ $\langle \alpha $j$|\beta $j$\rangle $ – $\langle \alpha $j$|$j$\beta \rangle $}\]

则需要知道解的信息。我们进一步将 \[VerticalSeparator]j\[RightAngleBracket] 展开，

\[\tilde{G}_{\alpha \beta }\text{ = }\sum _{j=1}^N \text{ $\langle \alpha $j$|\beta $j$\rangle $ – $\langle \alpha $j$|$j$\beta \rangle $ = }\sum _{j=1}^N \sum _{\gamma ,\delta =1}^K \text{ $\langle \gamma |$j$\rangle \langle $j$|\delta \rangle $($\langle \alpha \gamma |\beta \delta \rangle $ – $\langle \alpha \gamma |\delta \beta \rangle $) = }\sum _{\gamma ,\delta =1}^K \sum _{j=1}^N \text{ $\langle \gamma |$j$\rangle \langle $j$|\delta \rangle $($\langle \alpha \gamma |\beta \delta \rangle $ – $\langle \alpha \gamma |\delta \beta \rangle $) = }\sum _{\gamma ,\delta =1}^K \text{ (}\sum _{j=1}^N \text{ $\langle \gamma |$j$\rangle \langle $j$|\delta \rangle $) ($\langle \alpha \gamma |\beta \delta \rangle $ – $\langle \alpha \gamma |\delta \beta \rangle $) = }\sum _{\gamma ,\delta =1}^K D_{\gamma \delta }\text{($\langle \alpha \gamma |\beta \delta \rangle $ – $\langle \alpha \gamma |\delta \beta \rangle $) = }\sum _{\gamma ,\delta =1}^K G_{\alpha \beta \gamma \delta }D_{\gamma \delta }\text{ = (G$\cdot $D})_{\alpha \beta }\]

G\[CenterDot]D 是四维张量 G 和矩阵 D 的缩并 (contraction)。

\[G_{\alpha \beta \gamma \delta }\text{ = $\langle \alpha \gamma |\beta \delta \rangle $ – $\langle \alpha \gamma |\delta \beta \rangle $}\]

\[D_{\gamma \delta }\text{ = }\sum _{j=1}^N \text{ $\langle \gamma |$j$\rangle \langle $j$|\delta \rangle $ = }\sum _{j=1}^N C_{\text{$\gamma $j}}C_{\text{j$\delta $}}\text{ = }\sum _{j=1}^N C_{\text{$\gamma $j}}C_{\text{$\delta $j}}^*\]

注意 N < K，故 Subscript[D, \[Gamma]\[Delta]] 无法写成 C 中列向量的内积的形式，不过，如果我们对 C 进行截断，即只取其前 N 行构成 N*K 矩阵 B, 则有

\[D_{\gamma \delta }\text{ = }B_{\gamma }^{\dagger }B_{\delta }\]

写成这种形式结果上没有什么区别，但是尽量将求和写为矩阵或向量乘法对于实际的计算机处理非常有帮助，因为这意味着更加容易并行化或者更加容易利用先进的指令集。

一般来讲我们在量子力学中习惯将波函数展开系数列向量的内积构成的矩阵称为密度矩阵，因为它对应着密度算符在该基组下的矩阵表示

\[\hat{\rho }\text{ = $|\psi \rangle \langle \psi |$}\]

\[\langle \alpha |\hat{\rho }\text{$|\beta \rangle $ = $\langle \alpha |\psi \rangle \langle \psi |\beta \rangle $ = }c_{\alpha }c_{\beta }^*\]

现在我们可以看出，Fock矩阵可以写成

\[\text{F = }H^{\text{core}}\text{ + }\tilde{G}\text{ = T + }V^{\text{ne}}\text{ + G$\cdot $D}\]

我们注意到，以上的项中，只有 D 依赖于解，或者说依赖于组合系数矩阵 C。如果我们有了一个 C，则根据上式我们就能计算出 F，再用这个 F，我们又能解出新的 C。反复这个过程，直到 C 计算出的 F 解出的 C’ 与 C 相等为止。此时，Fock 算符中的等效势能项，刚好与由所有波函数生成的势场相同，这就是所谓的自洽条件，因此，我们这个方法被叫做自洽场方法。

以上的自洽场方法在数学上有很多非常严重的问题。比方说，其收敛是绝对没有任何保证的。实际上我们也没有说它一定能收敛。恰恰相反，以上的过程不仅经常难以收敛，有时候甚至会不停震荡或者干脆跑飞。再比方说，我们实际上不仅不能保证它收敛，而且连判断它能否收敛也是有困难的。也就是说，在迭代过程中，我们几乎无法知道继续计算下去到底能不能得到结果。

不过，为了不跑题，我们现在先假装不存在这个问题。等以后我们再来看怎样改善这个非常麻烦的情况。(可以先参考Jensen书101页大致了解一下几种方法的思路)

解 FC = SCE 在基组正交的条件下会更加容易。在基组正交的条件下，Overscript[S, ~] = 1，方程化为 Overscript[F, ~]Overscript[C, ~] = Overscript[C, ~]E。我们现在再来看如何将方程化为该形式。设

\[\text{C = X}\tilde{C}\]

注意 X 可能不是酉的。现在，我们有

\[\text{FC = SCE $\Rightarrow $ FX}\tilde{C}\text{ = SX}\tilde{C}\text{E $\Rightarrow $ }X^{\dagger }\text{FX}\tilde{C}\text{ = }X^{\dagger }\text{SX}\tilde{C}E\]

比较上式与 Overscript[F, ~]Overscript[C, ~] = Overscript[C, ~]E ，我们有

\[\text{1 = }X^{\dagger }\text{SX}\]

\[\tilde{F}\text{ = }X^{\dagger }\text{FX}\]

因此，我们可以随意选择任何使 1 = X^\[Dagger]SX 成立的矩阵 X 来正交化我们的方程。X 的选择有两种流行的方法：

1. 利用 S 是厄米的，因此 S^(1/2) 也是厄米的，故 S^(1/2)^\[Dagger] = S^(-1/2)。故取 X = S^(1/2) 即可。这被称为对称正交化 (symmetric orthogonalization)

2. 对角化 S，得到 S = U^\[Dagger]sU。取 X = Us^(1/2)，则

\[X^{\dagger }\text{SX = }s^{-1/2}U^{\dagger }\text{SUs}^{1/2}\text{ = }s^{-1/2}\text{ss}^{1/2}\text{ = 1}\]

这被称为正则正交化 (canonical orthogonalization)

以上两种方法理论上来说并没有任何区别，因为计算 S^(1/2) 一般来说要依赖于计算 s^(1/2) 。不过，如果基组的重叠非常大，在实际的计算中，计算 s^(-1/2) 可能会发生除以一个极小的数的问题，甚至在极端情况下，如果基组中存在线性相关的基组，则S的本征值中将会有一些是0，从而计算 s^(-1/2) 将会发生 Division By Zero 的错误。

如果我们预期会发生以上的情况，我们应该采用第二种方法。虽然看起来第二种方法并没有解决以上的问题，不过，在计算出s后，由于本征值的排列是随意的，我们大可以重新排列这些本征值，并且将可能造成数值计算问题的项进行截断处理。这实际上等价于抛弃那些线性相关或者近似线性相关的基组(反正包含它们对计算精度提升的贡献也并不很大)。

另外，我们也可以选择施密特正交化等大家耳熟能详的方法，不过，Schmidt 正交化的矩阵形式比较繁琐，而且由于不断投影，倾向于生成过大或者过小的矢量，通过在正交化的过程中不断重新归一可以解决这个问题。其它的正交化方法，可以参考高等代数。

到这里，我们已经大致讲完了 SCF 方法的整个推导过程。现在，我们把整个 SCF 的过程列出来总结一下，为我们接下来写程序提供一个思路。

1. 首先，我们选定一个要计算的分子。根据 BO 近似，我们固定核的坐标在某些位置，来计算电子能量。

2. 我们需要知道体系有几个电子。我们用 N 来代表体系的总电子数。

3. 我们选定一个基组。基组是任意的一系列电子自旋轨道。理论上来说，选择任何有3个空间自由度和1个自旋自由度的波函数都是允许的。我们用 K 来代表我们选择的基组的大小。

4. 计算 T, V^ne, G，X。这些在之后都不会再改变了。

5. 我们任意给定一个初始的密度猜测 C \[LeftArrow] Subscript[C, 0]

6. 备份现在的密度 C， Subscript[C, last] \[LeftArrow] C

7. 用 C 计算 D

8. 用 D 和 T, V^ne, G 计算 Fock 矩阵 F

\[\text{9. }\tilde{F}\text{ $\leftarrow $ }X^{\dagger }\text{FX}\]

10. 对角化 Overscript[F, ~] 得到 Overscript[F, ~]Overscript[C, ~] = Overscript[C, ~]E

\[\text{11. C $\leftarrow $ X}\tilde{C}\]

12. C 与 Subscript[C, last] 相比较是否相同？如果不同，回到第6步。

现在我们终于打算来写这个程序了。不过在写程序之前，我们还有最后一个问题没有解决，即到底如何计算 T, V^ne, G，X。这属于所谓的量子化学积分计算问题，我们现在来逐一解决。由于我们开始时真正定义的动量算符、势能算符和电子排斥等等实际上都只包含空间坐标的信息，故我们需要将轨道的自旋部分先处理掉，再来继续积分。

我们首先来考虑 restricted closed-shell 体系。 这一条件的意思是体系有 N 个占据轨道，N为偶数，这N个轨道是由 N/2 个空间轨道乘以相反的自旋函数得到。我们用 \[Phi](r) 表示空间轨道。

在 restricted closed-shell 条件下，将自旋轨道的自旋部分相互配对积分出去，则Fock算符简化为(此时库伦算符和交换算符的定义也改变为空间轨道形式)

\[\hat{F}_i\text{ = }\hat{h}_i\text{ + }\sum _{j=1}^{N/2} \text{ ( 2}\hat{J}_j\text{ + }\hat{K}_j\text{ )}\]

Hartree-Fock 方程简化为 Restricted Hartree-Fock (RHF) equation

\[\hat{F}_i\left|\phi _i\text{$\rangle $ = }\varepsilon _i\right|\phi _i\rangle\]

首先 T, V^ne, X 的计算没有难度。

\[T_{\alpha \beta }\text{ = $\langle \alpha |$-}\frac{1}{2}\nabla _i{}^2\text{$|\beta \rangle $ = $\langle $}\phi _{\alpha }\left(r_i\text{)$|$-}\frac{1}{2}\nabla _i{}^2\right|\phi _{\beta }(r_i\text{)$\rangle $ = -}\frac{1}{2}\int \phi _{\alpha }{}^*(r_i)\nabla _i{}^2\phi _{\beta }(r_i\text{)d}r_i\]

\[V_{\alpha \beta }^{\text{ne}}\text{ = $\langle \alpha |$}\sum _{a=1}^{N_{\text{nuc}}} \frac{1}{\left\left| R_a-r_i\right\right| }\text{$|\beta \rangle $ = $\langle $}\phi _{\alpha }\left(r_i\text{)$|$}\sum _{a=1}^{N_{\text{nuc}}} \frac{1}{\left\left| R_a-r_i\right\right| }\right|\phi _{\beta }(r_i\text{)$\rangle $ = $\int $}\phi _{\alpha }{}^*(r_i)\sum _{a=1}^{N_{\text{nuc}}} \frac{1}{\left\left| R_a-r_i\right\right| }\phi _{\beta }(r_i\text{)d}r_i\]

\[S_{\alpha \beta }\text{ = $\langle \alpha |\beta \rangle $ = $\langle $}\phi _{\alpha }\left(r_i\text{)$|$}\phi _{\beta }\right(r_i\text{)$\rangle $ = $\int $}\phi _{\alpha }{}^*(r_i)\phi _{\beta }(r_i\text{)d}r_i\]

\[\text{X = }S^{1/2}\]

以上的计算，一共是 K^2 个积分的复杂度

G 是四维张量，由于 \[LeftAngleBracket]\[Alpha]\[Beta]\[VerticalSeparator]\[Gamma]\[Delta]\[RightAngleBracket] 在 RHF 下具有八重对称性，因此我们只要计算 K^4/8 个积分即可

\[G_{\alpha \beta \gamma \delta }\text{ = $\langle \alpha \gamma |\beta \delta \rangle $ – $\langle \alpha \gamma |\delta \beta \rangle $ = $\int \int $}\phi _{\alpha }\left(r_i\right)\phi _{\gamma }(r_j)\hat{g}_{\text{ij}}(r_i\text{, }r_j)\phi _{\beta }(r_i)\phi _{\delta }(r_j\text{)d}r_idr_j\text{ – $\int \int $}\phi _{\alpha }(r_i)\phi _{\gamma }(r_j)\hat{g}_{\text{ij}}(r_i\text{, }r_j)\phi _{\delta }(r_i)\phi _{\beta }(r_j\text{)d}r_idr_j\text{ = $\int \int $}\phi _{\alpha }(r_i)\phi _{\gamma }(r_j)\frac{1}{\left\left| r_i-r_j\right\right| }\phi _{\beta }(r_i)\phi _{\delta }(r_j\text{)d}r_idr_j\text{ – $\int \int $}\phi _{\alpha }(r_i)\phi _{\gamma }(r_j)\frac{1}{\left\left| r_i-r_j\right\right| }\phi _{\delta }(r_i)\phi _{\beta }(r_j\text{)d}r_idr_j\text{ = $\int \int $}\phi _{\alpha }(r_i)\phi _{\beta }(r_i)\frac{1}{\left\left| r_i-r_j\right\right| }\phi _{\gamma }(r_j)\phi _{\delta }(r_j\text{)d}r_idr_j\text{ – $\int \int $}\phi _{\alpha }(r_i)\phi _{\delta }(r_i)\frac{1}{\left\left| r_i-r_j\right\right| }\phi _{\gamma }(r_j)\phi _{\beta }(r_j\text{)d}r_idr_j\text{ = ($\alpha \beta |\gamma \delta $) – ($\alpha \delta |\gamma \beta $)}\]

注意

\[\text{(12$|$34) = (12$|$43) = (21$|$43) = (21$|$34) = (34$|$21) = (34$|$12) = (43$|$12) = (43$|$21)}\]

有了以上的积分，我们就可以方便地进行 SCF 迭代了。方程为

\[\text{FC = SCE}\]

此方程称为 Roothaan 方程。在 unrestricted open-shell 体系下，则得到 Pople-Nesbet 方程。

现在我们试着按我们之前说的流程，写一个简单的程序看看。

电子耦合方法

基组

真正动手

一个简单的程序

Intel Math Kernel Library

优化

密度泛函理论(Density Functional Theory)

价键理论(Valence Bond Theory)

相对论性(Relativity)