Vasily's Blog

一个记录学习经历的站点

0%

使用SCF计算He原子基态能量

本文通过Hartree-Fock自洽场方法求出了He原子的基态能量。阅读本文需要读者对于Hartree-Fock自洽场方法有最基本的了解。在进行正式计算之前,我们首先明确在特定分子构型下,使用自洽场计算的过程。

如果读者觉得其中的某些积分计算过于困难,不妨使用Mathematica辅助积分。

一定分子构型下的自洽场计算过程

在分子结构一定时,自洽场分子轨道的计算过程如下:

  1. 选择一组基函数\(\chi_{s}\)

  2. 预估\(H_{rs}^{\mathrm{core}}\)\(S_{rs}\)和双电子积分\((rs|tu)\)

  3. 使用重叠积分\(S_{rs}\)和正交化计算矩阵\(A\)的元素,得到正交基函数

  4. 初猜\(\phi = \sum\limits_{s}c_{si}\chi_{s}\)\(c_{ci}\),并计算密度矩阵\(P\)

  5. 使用\(F_{rs} = H_{rs}^{\mathrm{core}} + \sum\limits_{t = 1}^{b}\sum\limits_{u = 1}^{b}P_{tu}[(rs|tu) - \dfrac{1}{2}(ru|ts)]\),并结合之前得到的密度矩阵\(P\),计算\(F_{rs}\)\((rs|tu)\)\(H_{rs}^{\mathrm{core}}\)

  6. 计算\(F' = A^{\dagger}FA\)

  7. 利用对角化计算\(F'\)对应的\(C'\)\(\varepsilon\)

  8. 计算\(C = AC'\)

  9. 计算新的密度矩阵\(P^{*} = 2CC^{\dagger}\)

  10. 检查新的密度矩阵\(P^{*}\)与原先的密度矩阵\(P\)的差别。如果二者较为接近,那么停止计算;如果二者相差较大,那么返回第5步,使用\(P^{*}\)重新计算,直至计算完成。

使用SCF计算He原子能量

初始化波函数

有关基函数与基组的知识在这里暂且不再赘述。我们对每个原子轨道使用2个Slater型轨道(Slater Type Orbital, STO)进行模拟。对于He原子而言,我们暂且只计算基态。STO的通式为 \[ \dfrac{1}{\sqrt{2n!}}\left( \dfrac{2\zeta}{a_{0}}\right)^{\frac{2n+1}{2}}r^{n-1}\mathrm{e}^{-\frac{\zeta r}{a_{0}}}Y_{l}^{m}(\phi, \theta) \] 对于基态而言,\(n = 1\)\(l = m = 0\)。将其代入,并使用原子单位可得 \[ \dfrac{1}{\sqrt{2}}(2\zeta)^{\frac{3}{2}}\mathrm{e}^{-\zeta r}Y_{0}^{0} \]\(Y_{0}^{0} = \dfrac{1}{\sqrt{4\pi}}\)代入可得 \[ \sqrt{\dfrac{\zeta^{3}}{\pi}}\mathrm{e}^{-\zeta r} \] 由此我们可以得到两个STO \[ \left \{ \begin{aligned} & \chi_{1} = \sqrt{\dfrac{\zeta_{1}^{3}}{\pi}}\mathrm{e}^{-\zeta_{1} r} \\ & \chi_{2} = \sqrt{\dfrac{\zeta_{2}^{3}}{\pi}}\mathrm{e}^{-\zeta_{2} r} \\ \end{aligned} \right. \] 此时假设\(\zeta_{1} = 1.45\)\(\zeta_{2} = 2.91\)。 ## 计算\(H_{rs}^{\mathrm{core}}\)\(S_{rs}\)以及\((rs|tu)\) 接下来计算\(H_{rs}^{\mathrm{core}}\)\(S_{rs}\)以及双电子积分\((rs|tu)\)。首先计算重叠积分\(S_{rs}\)。我们此时有3个需要计算的积分:\(S_{11}\)\(S_{12}\)\(S_{22}\)\(S_{11}\)\(S_{22}\)本质上相同,我们以\(S_{11}\)为例计算: \[ \begin{aligned} S_{11} = \langle \chi_{1} | \chi_{1} \rangle & = \dfrac{\zeta_{1}^{3}}{\pi}\iiint \mathrm{e}^{-2\zeta_{1} r} r^{2}\sin\phi\mathrm{d} r\mathrm{d} \phi \mathrm{d} \theta\\ & = \dfrac{\zeta_{1}^{3}}{\pi}\int_{0}^{2\pi} \mathrm{d} \theta \int_{0}^{\pi} \sin\phi \mathrm{d} \phi \int_{0}^{\infty} r^{2}\mathrm{e}^{-2\zeta_{1} r}\mathrm{d} r \\ & = 1 \end{aligned} \]

\(S_{12}\)等于 \[ \begin{aligned} S_{12} = \langle \chi_{1}|\chi_{2}\rangle & = \dfrac{\zeta^{\frac{3}{2}}_{1}\zeta_{2}^{\frac{3}{2}}}{\pi}\iiint \mathrm{e}^{-(\zeta_{1}+\zeta_{2}) r} r^{2}\sin\phi\mathrm{d} r\mathrm{d} \phi \mathrm{d} \theta\\ & = \dfrac{\zeta_{1}^{3}}{\pi}\int_{0}^{2\pi} \mathrm{d} \theta \int_{0}^{\pi} \sin\phi \mathrm{d} \phi \int_{0}^{\infty} r^{2}\mathrm{e}^{-(\zeta_{1}+\zeta_{2}) r}\mathrm{d} r \\ & = \dfrac{8\zeta_{1}^{\frac{3}{2}}\zeta_{2}^{\frac{3}{2}}}{(\zeta_{1} + \zeta_{2})^{3}} = 0.8366 \end{aligned} \]

然后计算\(H_{rs}^{\mathrm{core}}\)。我们有4个需要计算的积分:\(H_{11}^{\mathrm{core}}\)\(H_{12}^{\mathrm{core}}\)\(H_{21}^{\mathrm{core}}\)\(H_{22}^{\mathrm{core}}\)。首先计算\(H_{11}^{\mathrm{core}}\) \[ \begin{aligned} H_{11}^{\mathrm{core}} = \langle \chi_{1} | \hat{H}^{\mathrm{core}} | \chi_{1}\rangle & = \dfrac{\zeta_{1}^{3}}{\pi}\iiint \mathrm{e}^{-\zeta_{1} r} \left( \dfrac{1}{2}\nabla^{2} - \dfrac{2}{r}\right) \mathrm{e}^{-\zeta_{1} r}r^{2}\sin\phi\mathrm{d} r\mathrm{d} \phi \mathrm{d} \theta \\ & = \dfrac{\zeta_{1}^{3}}{\pi}\iiint \mathrm{e}^{-\zeta_{1} r} \int_{0}^{2\pi} \mathrm{d} \theta \int_{0}^{\pi} \sin\phi \mathrm{d} \phi\int_{0}^{\infty}-\dfrac{r}{2}\mathrm{e}^{-2\zeta_{1}r}(4 - 2\zeta_{1} + r\zeta_{1}^{2}) \\ & = \dfrac{1}{2}\zeta_{1}^{2} - 2\zeta_{1} = -1.8488 \end{aligned} \] 同理,\(H_{22}^{\mathrm{core}} = \dfrac{1}{2}\zeta_{2}^{2} - 2\zeta_{2} = -1.5860\)。而对于\(H_{12}^{\mathrm{core}}\)\(H_{21}^{\mathrm{core}}\),由于二者满足\(H_{12}^{\mathrm{core}} = H_{21}^{\mathrm{core}*}\),且二者均为实数,因此\(H_{12}^{\mathrm{core}} = H_{21}^{\mathrm{core}}\)。我们以\(H_{12}^{\mathrm{core}}\)为例进行计算 \[ \begin{aligned} H_{12}^{\mathrm{core}} = \langle \chi_{1}|\hat{H}|\chi_{2}\rangle & = \dfrac{\zeta^{\frac{3}{2}}_{1}\zeta_{2}^{\frac{3}{2}}}{\pi}\iiint \mathrm{e}^{-\zeta_{1} r} \left( \dfrac{1}{2}\nabla^{2} - \dfrac{2}{r}\right) \mathrm{e}^{-\zeta_{2} r}r^{2}\sin\phi\mathrm{d} r\mathrm{d} \phi \mathrm{d} \theta \\ & = \dfrac{\zeta^{\frac{3}{2}}_{1}\zeta_{2}^{\frac{3}{2}}}{\pi}\int_{0}^{2\pi} \mathrm{d} \theta \int_{0}^{\pi} \sin\phi \mathrm{d} \phi\int_{0}^{\infty}-\dfrac{r}{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r}(4 - 2\zeta_{2} + r\zeta_{2}^{2})\mathrm{d} r \\ & = \dfrac{\zeta^{\frac{3}{2}}_{1}\zeta_{2}^{\frac{3}{2}}(4\zeta_{1}\zeta_{2} - 8\zeta_{1} - 8\zeta_{2})}{(\zeta_{1} + \zeta_{2})^{3}} = -1.8826 \end{aligned} \]

最后计算双电子积分,回忆一下,双电子积分的定义如下 \[ (rs|tu) = \iint\dfrac{\chi_{r}^{*}(1)\chi_{s}(1)\chi_{t}^{*}(2)\chi_{u}^{*}(2)}{r_{12}} \mathrm{d} v_{1}\mathrm{d} v_{2} \] 由此,双电子积分的组合如下(下面只是一种表示方式,并不是矩阵的意思) \[ \begin{bmatrix} (11|11) & (11|12) & (11|21) & (11|22) \\ (12|11) & (12|12) & (12|21) & (12|22) \\ (21|11) & (21|12) & (21|21) & (21|22) \\ (22|11) & (22|12) & (22|21) & (22|22) \\ \end{bmatrix} \] 由于我们此时使用的是实函数,因此 \[ \begin{aligned} & (11|12) = (11|21) = (12|11) = (21|11)\\ & (12|12) = (12|21) = (21|12) = (21|21)\\ & (12|22) = (21|22) = (22|12) = (22|21)\\ & (11|22) = (22|11)\\ \end{aligned} \] 我们来一个个进行计算。首先计算\((11|11)\) \[ \begin{aligned} (11|11) & = \dfrac{\zeta_{1}^{6}}{\pi^{2}}\iint \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{1}r_{2}}}{r_{12}}\mathrm{d} v_{1}\mathrm{d} v_{2}\\ & = \dfrac{\zeta_{1}^{6}}{\pi^{2}}\int_{0}^{2\pi}\mathrm{d} \theta_{1} \int_{0}^{2\pi} \mathrm{d} \theta_{2} \int_{0}^{\pi} \sin\phi_{1} \mathrm{d} \phi_{1}\int_{0}^{\pi}\sin \phi_{2}\mathrm{d} \phi_{2}\int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{1}r_{2}}}{r_{12}}r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2} \\ & = 16\zeta_{1}^{6} \int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{1}r_{2}}}{r_{12}}r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2} \\ \end{aligned} \] 下面计算\(\displaystyle\int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{1}r_{2}}}{r_{12}}r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2}\),不难看出,重点是如何将\(r_{12}\)\(r_{1}\)\(r_{2}\)表示。我们在此直接给出结论,\(r_{12}\)可以表示为 \[ \dfrac{1}{r_{12}} = \sum\limits_{l = 0}^{\infty}\sum\limits_{m = -l}^{l}\dfrac{4\pi}{2l + 1}\dfrac{r_{\mathrm{min}}^{l}}{r_{\mathrm{max}}^{l + 1}}[Y_{l}^{m}(\phi_{1}, \theta_{1})]^{*}Y_{l}^{m}(\phi_{2}, \theta_{2}) \] 我们先不解释\(r_{\mathrm{max}}\)\(r_{\mathrm{min}}\)的含义。将\(l = m = 0\)代入可得 \[ \dfrac{1}{r_{12}} = \dfrac{1}{r_{\mathrm{max}}} \] 由此 \[ \int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{1}r_{2}}}{r_{\mathrm{max}}}r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2} \] 我们可以将\(r_{1}\)的积分区间分为两部分:\([0, r_{2}]\)\([r_{2}, \infty)\)。此时\(r_{\mathrm{max}}\)的含义为 \[ r_{\mathrm{max}} = \left \{ \begin{aligned} & r_{2},\quad r_{1} \in [0, r_{2}] \\ & r_{1},\quad r_{1} \in [r_{1}, \infty) \\ \end{aligned} \right. \] 由此 \[ \begin{aligned} & \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-2\zeta_{1}r_{2}}\mathrm{d} r_{2} \left( \int_{0}^{r_{2}}\dfrac{r_{1}^{2}\mathrm{e}^{-2\zeta_{1}r_{1}}}{r_{2}}\mathrm{d} r_{1} + \int_{r_{2}}^{\infty}r_{1}\mathrm{e}^{-2\zeta_{1}r_{1}} \mathrm{d} r_{1}\right) \\ & = \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-2\zeta_{1}r_{2}}\mathrm{d} r_{2} \left( \int_{0}^{r_{2}}\dfrac{r_{1}^{2}\mathrm{e}^{-2\zeta_{1}r_{1}}}{r_{2}}\mathrm{d} r_{1} + \int_{r_{2}}^{\infty}r_{1}\mathrm{e}^{-2\zeta_{1}r_{1}} \mathrm{d} r_{1}\right) \text{(用分部积分法)} \\ & = \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-2\zeta_{1}r_{2}}\mathrm{d} r_{2} \left[ \dfrac{1 - \mathrm{e}^{-2\zeta_{1}r_{2}}(1 + 2\zeta_{1}r_{2} + 2\zeta_{1}^{2}r_{2}^{2})}{4\zeta_{1}^{3}r_{2}} + \dfrac{\mathrm{e}^{-2\zeta_{1}r_{2}}(1 + 2\zeta_{1}r_{2})}{4\zeta_{1}^{2}} \right] \\ & = \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-2\zeta_{1}r_{2}}}{4\zeta_{1}^{3}}\mathrm{d} r_{2} - \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-4\zeta_{1}r_{2}}(1 + 2\zeta_{1}r_{2} + 2\zeta_{1}r_{2}^{2})}{4\zeta_{1}^{3}} \mathrm{d} r_{2} + \int_{0}^{\infty}\dfrac{r_{2}^{2}\mathrm{e}^{-4\zeta_{1}r_{2}}(1 + 2\zeta_{1}r_{2})}{4\zeta_{1}^{2}} \mathrm{d} r_{2} \\ & = \dfrac{5}{128\zeta_{1}^{5}} \\ \end{aligned} \] 由此可得 \[ \begin{aligned} (11|11) & = 16\zeta_{1}^{6} \cdot \dfrac{5}{128\zeta_{1}^{5}} \\ & = \dfrac{5\zeta_{1}}{8} = 0.9062\\ \end{aligned} \] 同理,\((22|22) = \dfrac{5\zeta_{2}}{8} = 1.8188\)

下面计算\((11|22)\) \[ \begin{aligned} (11|22) & = \dfrac{\zeta_{1}^{3}\zeta_{2}^{3}}{\pi^{2}}\iint \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{2}r_{2}}}{r_{12}}\mathrm{d} v_{1}\mathrm{d} v_{2}\\ & = 16\zeta_{1}^{3}\zeta_{2}^{3}\int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{2}r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2} \mathrm{d} r_{1}\mathrm{d} r_{2} \end{aligned} \] 针对\(\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{2}r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2} \mathrm{d} r_{1}\mathrm{d} r_{2}\),我们有 \[ \begin{aligned} & \int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-2\zeta_{2}r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2} \mathrm{d} r_{1}\mathrm{d} r_{2} \\ & = \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-2\zeta_{2}r_{2}}}{4\zeta_{1}^{3}}\mathrm{d} r_{2} - \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-2(\zeta_{1} + \zeta_{2})r_{2}}(1 + 2\zeta_{1}r_{2} + 2\zeta_{1}r_{2}^{2})}{4\zeta_{1}^{3}} \mathrm{d} r_{2} + \int_{0}^{\infty}\dfrac{r_{2}^{2}\mathrm{e}^{-2(\zeta_{1} + \zeta_{2})r_{2}}(1 + 2\zeta_{1}r_{2})}{4\zeta_{1}^{2}} \mathrm{d} r_{2} \\ & = \dfrac{1}{16\zeta_{1}^{3}\zeta_{2}^{2}} - \dfrac{6\zeta_{1}^{2} + 4\zeta_{1}\zeta_{2} + \zeta_{2}^{2}}{16\zeta_{1}^{3}(\zeta_{1} + \zeta_{3})^{4}} + \dfrac{4\zeta_{1} + \zeta_{2}}{16\zeta_{1}^{2}(\zeta_{1} + \zeta_{2})^{4}} \\ & = \dfrac{\zeta_{1}^{4}\zeta_{2} + 4\zeta_{1}^{3}\zeta_{2}^{2} + 4\zeta_{1}^{2}\zeta_{2}^{3} + \zeta_{1}\zeta_{2}^{4}}{16\zeta_{1}^{3}\zeta_{2}^{3}(\zeta_{1} + \zeta_{2})^{4}} \\ \end{aligned} \] 由此 \[ \begin{aligned} (11|22) & = 16\zeta_{1}^{3}\zeta_{2}^{3} \cdot \dfrac{\zeta_{1}^{4}\zeta_{2} + 4\zeta_{1}^{3}\zeta_{2}^{2} + 4\zeta_{1}^{2}\zeta_{2}^{3} + \zeta_{1}\zeta_{2}^{4}}{16\zeta_{1}^{3}\zeta_{2}^{3}(\zeta_{1} + \zeta_{2})^{4}} \\ & = \dfrac{\zeta_{1}^{4}\zeta_{2} + 4\zeta_{1}^{3}\zeta_{2}^{2} + 4\zeta_{1}^{2}\zeta_{2}^{3} + \zeta_{1}\zeta_{2}^{4}}{(\zeta_{1} + \zeta_{2})^{4}} = 1.1826 \\ \end{aligned} \]

下面计算\((11|12)\) \[ \begin{aligned} (11|12) & = \dfrac{\zeta_{1}^{\frac{9}{2}}\zeta_{2}^{\frac{3}{2}}}{\pi^{2}}\iint \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}}\mathrm{d} v_{1}\mathrm{d} v_{2}\\ & = 16\zeta_{1}^{\frac{9}{2}}\zeta_{2}^{\frac{3}{2}} \int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2} \mathrm{d} r_{1} \mathrm{d} r_{2} \\ \end{aligned} \] 针对\(\displaystyle \int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2} \mathrm{d} r_{1} \mathrm{d} r_{2}\),我们有 \[ \begin{aligned} & \int_{0}^{\infty}\int_{0}^{\infty} \dfrac{\mathrm{e}^{-2\zeta_{1}r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2} \mathrm{d} r_{1} \mathrm{d} r_{2} \\ & = \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}\mathrm{d} r_{2} \left( \int_{0}^{r_{2}}\dfrac{r_{1}^{2}\mathrm{e}^{-2\zeta_{1}r_{1}}}{r_{2}}\mathrm{d} r_{1} + \int_{r_{2}}^{\infty}r_{1}\mathrm{e}^{-2\zeta_{1}r_{1}} \mathrm{d} r_{1}\right) \\ & = \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}\mathrm{d} r_{2} \left[ \dfrac{1 - \mathrm{e}^{-2\zeta_{1}r_{2}}(1 + 2\zeta_{1}r_{2} + 2\zeta_{1}^{2}r_{2}^{2})}{4\zeta_{1}^{3}r_{2}} + \dfrac{\mathrm{e}^{-2\zeta_{1}r_{2}}(1 + 2\zeta_{1}r_{2})}{4\zeta_{1}^{2}} \right] \\ & = \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{4\zeta_{1}^{3}}\mathrm{d} r_{2} - \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-(3\zeta_{1}+ \zeta_{2})r_{2}}(1 + 2\zeta_{1}r_{2} + 2\zeta_{1}^{2}r_{2}^{2})}{4\zeta_{1}^{3}} \mathrm{d} r_{2} + \int_{0}^{\infty}\dfrac{r_{2}^{2}\mathrm{e}^{-(3\zeta_{1} + \zeta_{2})r_{2}}(1 + 2\zeta_{1}r_{2})}{4\zeta_{1}^{2}} \mathrm{d} r_{2} \\ & = \dfrac{1}{4\zeta_{1}^{3}(\zeta_{1} + \zeta_{2})^{2}} - \dfrac{33\zeta_{1}^{2} + 10\zeta_{1}\zeta_{2} + \zeta_{2}^{2}}{4\zeta_{1}^{3}(3\zeta_{1} + \zeta_{2})^{4}} + \dfrac{9\zeta_{1} + \zeta_{2}}{2\zeta_{1}^{2}(3\zeta_{1} + \zeta_{2})^{4}}\\ & = \dfrac{1}{(3\zeta_{1} + \zeta_{2})^{4}}\left[ \dfrac{12\zeta_{1} + 8\zeta_{2}}{(\zeta_{1} + \zeta_{2})^{2}} + \dfrac{9\zeta_{1} + \zeta_{2}}{2\zeta_{1}^{2}} \right]\\ \end{aligned} \] 由此可得 \[ \begin{aligned} (11|12) & = 16\zeta_{1}^{\frac{9}{2}}\zeta_{2}^{\frac{3}{2}} \cdot \dfrac{1}{(3\zeta_{1} + \zeta_{2})^{4}}\left[ \dfrac{12\zeta_{1} + 8\zeta_{2}}{(\zeta_{1} + \zeta_{2})^{2}} + \dfrac{9\zeta_{1} + \zeta_{2}}{2\zeta_{1}^{2}} \right] \\ & = \dfrac{16\zeta_{1}^{\frac{9}{2}}\zeta_{2}^{\frac{3}{2}}}{(3\zeta_{1} + \zeta_{2})^{4}}\left[ \dfrac{12\zeta_{1} + 8\zeta_{2}}{(\zeta_{1} + \zeta_{2})^{2}} + \dfrac{9\zeta_{1} + \zeta_{2}}{2\zeta_{1}^{2}} \right] = 0.9033 \\ \end{aligned} \] 由此不难得到\((22|21)\)\[ (22|21) = \dfrac{\zeta_{2}^{\frac{9}{2}}\zeta_{1}^{\frac{3}{2}}}{(3\zeta_{2} + \zeta_{1})^{4}}\left[ \dfrac{12\zeta_{2} + 8\zeta_{1}}{(\zeta_{2} + \zeta_{1})^{2}} + \dfrac{9\zeta_{2} + \zeta_{1}}{2\zeta_{2}^{2}} \right] = 1.2980 \]

下面计算\((12|12)\) \[ \begin{aligned} (12|12) & = \dfrac{\zeta_{1}^{3}\zeta^{3}}{\pi^{2}} \iint \dfrac{\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}}\mathrm{d} v_{1}\mathrm{d} v_{2}\\ & = 16\zeta_{1}^{3}\zeta^{3}\int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2} \end{aligned} \] 针对\(\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2}\),我们有 \[ \begin{aligned} & \int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{1}}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{r_{12}} r_{1}^{2}r_{2}^{2}\mathrm{d} r_{1}\mathrm{d} r_{2} \\ & = \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}\mathrm{d} r_{2} \left( \int_{0}^{r_{2}}\dfrac{r_{1}^{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{1}}}{r_{2}}\mathrm{d} r_{1} + \int_{r_{2}}^{\infty}r_{1}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{1}} \mathrm{d} r_{1}\right) \\ & = \int_{0}^{\infty}r_{2}^{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}\mathrm{d} r_{2} \left[ \dfrac{2 - \mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}(2+2r_{2}(\zeta_{1} + \zeta_{2})+r_{2}^{2}(\zeta_{1} + \zeta_{2})^{2})}{r_{2}(\zeta_{1} + \zeta_{2})^{3}} + \dfrac{\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}(1 + r_{2}(\zeta_{1} + \zeta_{2}))}{(\zeta_{1} + \zeta_{2})^{2}} \right] \\ & = \int_{0}^{\infty}\dfrac{2r_{2}\mathrm{e}^{-(\zeta_{1} + \zeta_{2})r_{2}}}{(\zeta_{1} + \zeta_{2})^{3}}\mathrm{d} r_{2} - \int_{0}^{\infty}\dfrac{r_{2}\mathrm{e}^{-2(\zeta_{1}+ \zeta_{2})r_{2}}(2+2r_{2}(\zeta_{1} + \zeta_{2})+r_{2}^{2}(\zeta_{1} + \zeta_{2})^{2})}{(\zeta_{1} + \zeta_{2})^{3}} \mathrm{d} r_{2} \\ & \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad + \int_{0}^{\infty}\dfrac{r_{2}^{2}\mathrm{e}^{-2(\zeta_{1} + \zeta_{2})r_{2}}(1 + r_{2}(\zeta_{1} + \zeta_{2}))}{(\zeta_{1} + \zeta_{2})^{2}} \mathrm{d} r_{2} \\ & = \dfrac{2}{(\zeta_{1} + \zeta_{2})^{5}} - \dfrac{11}{8(\zeta_{1} + \zeta_{2})^{5}} + \dfrac{5}{8(\zeta_{1} + \zeta_{2})^{5}}\\ & = \dfrac{5}{4(\zeta_{1} + \zeta_{2})^{5}}\\ \end{aligned} \] 由此可得 \[ \begin{aligned} (12|12) & = 16\zeta_{1}^{3}\zeta_{2}^{3} \cdot \dfrac{5}{4(\zeta_{1} + \zeta_{2})^{5}} \\ & = \dfrac{20\zeta_{1}^{3}\zeta_{2}^{3}}{(\zeta_{1} + \zeta_{2})^{5}} = 0.9536 \\ \end{aligned} \]

总结一下,我们现在已经得到了如下的积分结果:

  1. \(S_{rs}\) \[ \left \{ \begin{aligned} & S_{11} = S_{22} = 1 \\ & S_{12} = S_{21} = 0.8366 \\ \end{aligned} \right. \]

  2. \(H_{rs}^{\mathrm{core}}\) \[ \left \{ \begin{aligned} & H_{11}^{\mathrm{core}} = -1.8488 \\ & H_{12}^{\mathrm{core}} = H_{21}^{\mathrm{core}} = -1.8826 \\ & H_{22}^{\mathrm{core}} = -1.5860\\ \end{aligned} \right. \]

  3. \((rs|tu)\) \[ \left \{ \begin{aligned} & (11|11) = 0.9062 \\ & (22|22) = 1.8188 \\ & (11|22) = (22|11) = 1.1826 \\ & (11|12) = (11|21) = (12|11) = (21|11) = 0.9033\\ & (12|22) = (21|22) = (22|12) = (22|21) = 1.2980\\ & (12|12) = (12|21) = (21|12) = (21|21) = 0.9536 \\ \end{aligned} \right. \]

给定初猜波函数

由于我们此时有两个基函数,因此初猜波函数应当为此二者的线性组合 \[ \phi_{1} = c_{11}\chi_{1} + c_{21}\chi_{2} \]

接下来我们要对\(c_{11}\)\(c_{21}\)进行赋值。由于波函数要满足归一化,因此重点不是两个系数的绝对数值大小,而是相对大小。对于Slater型基函数 \[ \chi = \sqrt{\dfrac{\zeta^{3}}{\pi}} \mathrm{e}^{-\zeta r} \] 其是从类氢原子体系的波函数演变而来的,因此其中\(\zeta\)表示类氢原子的有效核电荷数。显然,由于电子之间存在屏蔽效应,因此对于其中一个电子而言,其感受到的有效核电荷小于2。这样,\(\zeta_{1}\)对应的波函数\(\chi_{1}\)应当更加接近真实的波函数,从而要求其在线性组合中占有更大的比例。因此\(c_{11}\)应当比\(c_{21}\)大。我们暂且假定\(\dfrac{c_{11}}{c_{21}} = 2\)。由归一化得 \[ \begin{aligned} \int |\phi_{1}|^{2} \mathrm{d} \tau & = c_{11}^{2} + c_{21}^{2} + 2c_{11}c_{21}S_{12} \\ & = (5 + 4S_{12})c_{21}^{2} = 1\\ \end{aligned} \] 由此可得 \[ \left \{ \begin{aligned} & c_{11} = 0.6922 \\ & c_{21} = 0.3461 \\ \end{aligned} \right. \] 由此初猜波函数为 \[ \phi_{1} = 0.6922\chi_{1} + 0.3461\chi_{2} \]

计算密度矩阵

密度矩阵\(P\)中每个矩阵元的定义为 \[ P_{tu} = 2\sum\limits_{j = 1}^{\frac{n}{2}} c_{tj}^{*}c_{uj} \] 由此可得 \[ \left \{ \begin{aligned} & P_{11} = 2c_{11}^{2} = 0.9583 \\ & P_{12} = P_{21} = 2c_{11}c_{21} = 0.4791 \\ & P_{22} = 2c_{21}^{2} = 0.2396\\ \end{aligned} \right. \]

计算Fock矩阵元

根据\(F_{rs} = H_{rs}^{\mathrm{core}} + \sum\limits_{t = 1}^{b}\sum\limits_{u = 1}^{b}P_{tu}[(rs|tu) - \dfrac{1}{2}(ru|ts)]\),我们可以得到Fock矩阵的矩阵元 \[ \left \{ \begin{aligned} & F_{11} = H_{11}^{\mathrm{core}} + P_{11}[(11|11) - \dfrac{1}{2}(11|11)] + P_{12}[(11|12) - \dfrac{1}{2}(12|11)] \\ & \quad\quad\quad + P_{21}[(11|21) - \dfrac{1}{2}(11|21)] + P_{22}[(11|22) - \dfrac{1}{2}(12|21)] = -0.813\\ & F_{12} = F_{21} = H_{12}^{\mathrm{core}} + P_{11}[(12|11) - \dfrac{1}{2}(11|12)] + P_{12}[(12|12) - \dfrac{1}{2}(12|12)] \\ & \quad\quad\quad + P_{21}[(12|21) - \dfrac{1}{2}(11|22)] + P_{22}[(12|22) - \dfrac{1}{2}(12|22)] = -0.892\\ & F_{22} = H_{22}^{\mathrm{core}} + P_{11}[(22|11) - \dfrac{1}{2}(21|12)] + P_{12}[(22|12) - \dfrac{1}{2}(22|12)] \\ & \quad\quad\quad + P_{21}[(22|21) - \dfrac{1}{2}(21|22)] + P_{22}[(22|22) - \dfrac{1}{2}(22|22)] = -0.070\\ \end{aligned} \right. \]

计算系数矩阵\(C\)和轨道能矩阵\(\varepsilon\)

我们现在得到了给定基函数下的Fock矩阵 \[ F = \begin{bmatrix} -0.813 & -0.892 \\ -0.892 & -0.070 \\ \end{bmatrix} \] 由于Fock矩阵与系数矩阵、轨道能矩阵满足 \[ FC = SC\varepsilon \] 其中,\(S\)为矩阵元为重叠积分的矩阵,即 \[ S = \begin{bmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \\ \end{bmatrix} \]

由于我们现在使用的并不是正交基,因此要求出\(C\)\(\varepsilon\)有两种办法。我们先介绍一种较为通用(但此处我们不用)的方法,然后介绍一种对于小体系比较方便的方法。

由于基函数之间线性无关,因此我们总能用其构造出一组正交基。也就是说,存在一组基函数\(\{\chi_{i}'\}\)满足\(\chi_{i}' = \sum\limits_{j}a_{ji}\chi_{j}\),并且此时新的重叠矩阵\(S'\)满足 \[ S_{rs}' = \langle \chi_{r}'|\chi_{s}'\rangle = \delta_{rs} \] 此时波函数\(\phi_{i}\)可以表示为 \[ \phi_{i} = \begin{bmatrix} c_{1i}' & c_{2i}' & \dots & c_{ni}'\\ \end{bmatrix} \begin{bmatrix} \chi_{1}' \\ \chi_{2}' \\ \vdots \\ \chi_{n}' \\ \end{bmatrix} = \begin{bmatrix} c_{1i}' & c_{2i}' & \dots & c_{ni}'\\ \end{bmatrix} \begin{bmatrix} a_{11} & a_{21} & \dots & a_{n1} \\ a_{12} & a_{22} & \dots & a_{n2} \\ \vdots & \vdots & & \vdots \\ a_{1n} & a_{2n} & \dots & a_{nn} \\ \end{bmatrix} \begin{bmatrix} \chi_{1} \\ \chi_{2} \\ \vdots \\ \chi_{n} \\ \end{bmatrix} \] 我们将波函数由新的基函数线性表出的系数矩阵记为\(C'\),将原有基函数用新的基函数线性表出的系数矩阵记为\(A\),这样不难看出 \[ C^{\mathrm{T}} = C'^{\mathrm{T}}A^{\mathrm{T}} \] 从而 \[ C = AC' \]

此时我们将\(C' = CA\)代入\(FC = SC\varepsilon\)可得 \[ FAC' = SAC'\varepsilon \] 此时我们用\(A^{\dagger}\)来表示\((A^{*})^{\mathrm{T}}\)。两边同时左乘\(A^{\dagger}\)可得 \[ (A^{\dagger}FA)C' = (A^{\dagger}SA)C'\varepsilon \] 之所以前面那么别扭地假设\(\chi_{i}' = \sum\limits_{j}a_{ji}\chi_{j}\),而不是假设\(\chi_{i} = \sum\limits_{j}a_{ij}\chi_{j}'\)或者\(\chi_{i}' = \sum\limits_{j}a_{ij}\chi_{j}\),纯粹就是为了这里推导方便。如果矩阵\(A\)在矩阵\(C'\)或者\(C\)右侧,那么会得到类似 \(FC'A^{-1} = SC'A^{-1}\varepsilon\)的结果。此时\(A\)并没有直接和\(S\)接触,我们就无法进行进一步推导。

此时注意力集中的读者可能会想到,在新的基函数下Fock矩阵变为\(F' = A^{\dagger}FA\),而重叠矩阵变为\(S' = A^{\dagger}SA\)。由此可得 \[ F'C' = S'C'\varepsilon \] 但是我们不能只从物理上得出这一点,还要从数学上说明。对于\(S'\)的矩阵元\(S'_{rs}\),我们有(需要注意,新的基下必然存在\(S'\),我们这里是在寻找\(S'\)\(S\)的联系。也就是说,此时我们不知道\(S' = A^{\dagger}SA\)\[ \begin{aligned} S_{rs}' & = \langle \chi_{r}' | \chi_{s}'\rangle = \delta_{rs} \\ & = \sum\limits_{j}\sum\limits_{k} a_{jr}^{*}\langle \chi_{j} | \chi_{k}\rangle a_{ks} \\ & = \sum\limits_{j}\sum\limits_{k} a_{jr}^{*}S_{jk} a_{ks}\\ \end{aligned} \] 由矩阵元之间的关系不难得出 \[ A^{\dagger} S A = S' = I \]

接下来需要确定矩阵\(A\)。首先,由于矩阵\(S\)为实矩阵,并且是对称矩阵,因此\(S^{\dagger} = S\)。我们将满足这种性质的矩阵称为Hermite矩阵。在数学上可以证明,对于Hermite矩阵,总是存在酉矩阵\(U\),使得\(S\)可以被对角化 \[ U^{\dagger}SU = \Lambda \] 其中酉矩阵指的是满足\(U^{\dagger} = U^{-1}\)的矩阵。下面我们给出两种从\(U\)\(A\)的方法。

  1. 对称正交化:假设\(A = S^{-\frac{1}{2}}\)。为了构建这一矩阵,我们需要先将\(S\)对角化,然后将每个本征值倒数的平方根组成矩阵\(\Lambda^{-\frac{1}{2}}\),最后通过“逆对角化”得到\(S^{-\frac{1}{2}}\) \[ U\Lambda^{-\frac{1}{2}}U^{\dagger} = S^{-\frac{1}{2}} \] 不难证明\(ASA^{-1} = I\)。但当原来的基函数中有线性相关(或者接近线性相关)的情况时,这一方法在精度上会出现问题。因此我们更加推荐方法2。

  2. 正则正交化:此时我们假设\(A = U\Lambda^{-\frac{1}{2}}\)。由此 \[ A^{\dagger}SA = \Lambda^{-\frac{1}{2}}U^{\dagger}SU\Lambda^{-\frac{1}{2}} = I \] 这一方法在理论上精度更高。

无论如何,我们可以通过将矩阵\(S\)对角化的矩阵\(U\)来确定矩阵\(A\),从而由\(S' = A^{\dagger}SA\)确定\(S'\)。同时,我们也可以通过\(A\)确定\(F' = A^{\dagger}FA\)。由于\(S' = I\),因此原方程变为 \[ F'C' = C'\varepsilon \] 由于\(\varepsilon\)为对角矩阵,因此我们只需将\(F'\)对角化即可求出\(C'\)

总结上面的过程,我们可以得到

对于小体系,由于 \[ FC = SC\varepsilon \] 因此有 \[ F\boldsymbol{c}_{i} = \varepsilon_{i}S\boldsymbol{c}_{i} \] 假设\(F\)中每个行向量为\(\boldsymbol{f}^{\mathrm{T}}\)\(S\)中每个行向量为\(\boldsymbol{s}^{\mathrm{T}}\),那么有 \[ \boldsymbol{f}^{\mathrm{T}}\boldsymbol{c}_{i} = \varepsilon_{i}\boldsymbol{s}^{\mathrm{T}}\boldsymbol{c}_{i} \] 因此可以通过求\(\det(F - \varepsilon_{i}S) = 0\)来得到可能的所有\(\varepsilon_{i}\),从而由\((\boldsymbol{f}^{\mathrm{T}} - \varepsilon_{i}\boldsymbol{s}^{\mathrm{T}})\boldsymbol{c}_{i} = 0\)求出\(\boldsymbol{c}_{i}\),然后由\(\boldsymbol{c}_{i}\)组成矩阵\(C\)。由于我们这里只有一个波函数,因此使用这个方法。此时行列式为 \[ \begin{vmatrix} -0.813-\varepsilon_{i} & -0.892-0.8366\varepsilon_{i}\\ -0.892-0.8366\varepsilon_{i} & -0.07-\varepsilon_{i} \\ \end{vmatrix} = 0 \] 解方程可得 \[ \left \{ \begin{aligned} & \varepsilon_{1} = -0.853 \\ & \varepsilon_{2} = 2.884 \\ \end{aligned} \right. \] 我们选择能量较低的\(\varepsilon_{1}\)。由此可得 \[ c_{11}(F_{21} - \varepsilon_{1}S_{21}) + c_{21}(F_{22} - \varepsilon_{1}S_{22}) = 0 \] 由此可得 \[ \dfrac{c_{11}}{c_{21}} = 4.40 \] 由归一化不难得到 \[ c_{11}^{2} + c_{21}^{2} + 2c_{11}c_{21}S_{12} = 1 \]\(S_{12} = 0.8366\)以及\(\dfrac{c_{11}}{c_{21}} = 4.40\)代入可得 \[ \left \{ \begin{aligned} & c_{11} = 0.836 \\ & c_{21} = 0.190 \\ \end{aligned} \right. \]

第二次计算密度矩阵

此时,新密度矩阵的矩阵元分别为 \[ \left \{ \begin{aligned} & P_{11} = 2c_{11}^{2} = 1.398 \\ & P_{12} = P_{21} = 2c_{11}c_{21} = 0.318 \\ & P_{22} = 2c_{21}^{2} = 0.072 \\ \end{aligned} \right. \] 不难得出,此时密度矩阵与原有密度矩阵相差较大,因此我们需要重新计算Fock矩阵,进而重新求解\(C\)\(\varepsilon\)

第二次计算Fock矩阵元

根据\(F_{rs} = H_{rs}^{\mathrm{core}} + \sum\limits_{t = 1}^{b}\sum\limits_{u = 1}^{b}P_{tu}[(rs|tu) - \dfrac{1}{2}(ru|ts)]\),我们可以得到Fock矩阵的矩阵元 \[ \left \{ \begin{aligned} & F_{11} = H_{11}^{\mathrm{core}} + P_{11}[(11|11) - \dfrac{1}{2}(11|11)] + P_{12}[(11|12) - \dfrac{1}{2}(12|11)] \\ & \quad\quad\quad + P_{21}[(11|21) - \dfrac{1}{2}(11|21)] + P_{22}[(11|22) - \dfrac{1}{2}(12|21)] = -0.877\\ & F_{12} = F_{21} = H_{12}^{\mathrm{core}} + P_{11}[(12|11) - \dfrac{1}{2}(11|12)] + P_{12}[(12|12) - \dfrac{1}{2}(12|12)] \\ & \quad\quad\quad + P_{21}[(12|21) - \dfrac{1}{2}(11|22)] + P_{22}[(12|22) - \dfrac{1}{2}(12|22)] = -0.940\\ & F_{22} = H_{22}^{\mathrm{core}} + P_{11}[(22|11) - \dfrac{1}{2}(21|12)] + P_{12}[(22|12) - \dfrac{1}{2}(22|12)] \\ & \quad\quad\quad + P_{21}[(22|21) - \dfrac{1}{2}(21|22)] + P_{22}[(22|22) - \dfrac{1}{2}(22|22)] = -0.125\\ \end{aligned} \right. \]

第二次计算系数矩阵\(C\)和轨道能矩阵\(\varepsilon\)

求解行列式方程 \[ \begin{vmatrix} -0.877-\varepsilon_{i} & -0.940-0.8366\varepsilon_{i}\\ -0.940-0.8366\varepsilon_{i} & -0.125-\varepsilon_{i} \\ \end{vmatrix} = 0 \] 可得 \[ \left \{ \begin{aligned} & \varepsilon_{1} = -0.915 \\ & \varepsilon_{2} = 2.817 \\ \end{aligned} \right. \] 将能量较低的\(\varepsilon_{1}\)代入 \[ c_{11}(F_{21} - \varepsilon_{1}S_{21}) + c_{21}(F_{22} - \varepsilon_{1}S_{22}) = 0 \] 可得 \[ \dfrac{c_{11}}{c_{21}} = 4.53 \] 由归一化不难得到 \[ c_{11}^{2} + c_{21}^{2} + 2c_{11}c_{21}S_{12} = 1 \]\(S_{12} = 0.8366\)以及\(\dfrac{c_{11}}{c_{21}} = 4.53\)代入可得 \[ \left \{ \begin{aligned} & c_{11} = 0.840 \\ & c_{21} = 0.185 \\ \end{aligned} \right. \]

第三次计算密度矩阵

新密度矩阵的矩阵元分别为 \[ \left \{ \begin{aligned} & P_{11} = 2c_{11}^{2} = 1.411 \\ & P_{12} = P_{21} = 2c_{11}c_{21} = 0.311 \\ & P_{22} = 2c_{21}^{2} = 0.069 \\ \end{aligned} \right. \]

第三次计算Fock矩阵元

根据\(F_{rs} = H_{rs}^{\mathrm{core}} + \sum\limits_{t = 1}^{b}\sum\limits_{u = 1}^{b}P_{tu}[(rs|tu) - \dfrac{1}{2}(ru|ts)]\),我们可以得到Fock矩阵的矩阵元 \[ \left \{ \begin{aligned} & F_{11} = H_{11}^{\mathrm{core}} + P_{11}[(11|11) - \dfrac{1}{2}(11|11)] + P_{12}[(11|12) - \dfrac{1}{2}(12|11)] \\ & \quad\quad\quad + P_{21}[(11|21) - \dfrac{1}{2}(11|21)] + P_{22}[(11|22) - \dfrac{1}{2}(12|21)] = -0.880\\ & F_{12} = F_{21} = H_{12}^{\mathrm{core}} + P_{11}[(12|11) - \dfrac{1}{2}(11|12)] + P_{12}[(12|12) - \dfrac{1}{2}(12|12)] \\ & \quad\quad\quad + P_{21}[(12|21) - \dfrac{1}{2}(11|22)] + P_{22}[(12|22) - \dfrac{1}{2}(12|22)] = -0.940\\ & F_{22} = H_{22}^{\mathrm{core}} + P_{11}[(22|11) - \dfrac{1}{2}(21|12)] + P_{12}[(22|12) - \dfrac{1}{2}(22|12)] \\ & \quad\quad\quad + P_{21}[(22|21) - \dfrac{1}{2}(21|22)] + P_{22}[(22|22) - \dfrac{1}{2}(22|22)] = -0.124\\ \end{aligned} \right. \]

第三次计算系数矩阵\(C\)和轨道能矩阵\(\varepsilon\)

求解行列式方程 \[ \begin{vmatrix} -0.880-\varepsilon_{i} & -0.940-0.8366\varepsilon_{i}\\ -0.940-0.8366\varepsilon_{i} & -0.124-\varepsilon_{i} \\ \end{vmatrix} = 0 \] 可得 \[ \left \{ \begin{aligned} & \varepsilon_{1} = -0.917 \\ & \varepsilon_{2} = 2.813 \\ \end{aligned} \right. \] 将能量较低的\(\varepsilon_{1}\)代入 \[ c_{11}(F_{21} - \varepsilon_{1}S_{21}) + c_{21}(F_{22} - \varepsilon_{1}S_{22}) = 0 \] 可得 \[ \dfrac{c_{11}}{c_{21}} = 4.59 \] 由归一化不难得到 \[ c_{11}^{2} + c_{21}^{2} + 2c_{11}c_{21}S_{12} = 1 \]\(S_{12} = 0.8366\)以及\(\dfrac{c_{11}}{c_{21}} = 4.59\)代入可得 \[ \left \{ \begin{aligned} & c_{11} = 0.842 \\ & c_{21} = 0.183 \\ \end{aligned} \right. \]

第四次计算密度矩阵

新密度矩阵的矩阵元分别为 \[ \left \{ \begin{aligned} & P_{11} = 2c_{11}^{2} = 1.418 \\ & P_{12} = P_{21} = 2c_{11}c_{21} = 0.308 \\ & P_{22} = 2c_{21}^{2} = 0.067 \\ \end{aligned} \right. \]

第四次计算Fock矩阵元

根据\(F_{rs} = H_{rs}^{\mathrm{core}} + \sum\limits_{t = 1}^{b}\sum\limits_{u = 1}^{b}P_{tu}[(rs|tu) - \dfrac{1}{2}(ru|ts)]\),我们可以得到Fock矩阵的矩阵元 \[ \left \{ \begin{aligned} & F_{11} = H_{11}^{\mathrm{core}} + P_{11}[(11|11) - \dfrac{1}{2}(11|11)] + P_{12}[(11|12) - \dfrac{1}{2}(12|11)] \\ & \quad\quad\quad + P_{21}[(11|21) - \dfrac{1}{2}(11|21)] + P_{22}[(11|22) - \dfrac{1}{2}(12|21)] = -0.881\\ & F_{12} = F_{21} = H_{12}^{\mathrm{core}} + P_{11}[(12|11) - \dfrac{1}{2}(11|12)] + P_{12}[(12|12) - \dfrac{1}{2}(12|12)] \\ & \quad\quad\quad + P_{21}[(12|21) - \dfrac{1}{2}(11|22)] + P_{22}[(12|22) - \dfrac{1}{2}(12|22)] = -0.940\\ & F_{22} = H_{22}^{\mathrm{core}} + P_{11}[(22|11) - \dfrac{1}{2}(21|12)] + P_{12}[(22|12) - \dfrac{1}{2}(22|12)] \\ & \quad\quad\quad + P_{21}[(22|21) - \dfrac{1}{2}(21|22)] + P_{22}[(22|22) - \dfrac{1}{2}(22|22)] = -0.124\\ \end{aligned} \right. \]

第四次计算系数矩阵\(C\)和轨道能矩阵\(\varepsilon\)

求解行列式方程 \[ \begin{vmatrix} -0.881-\varepsilon_{i} & -0.940-0.8366\varepsilon_{i}\\ -0.940-0.8366\varepsilon_{i} & -0.124-\varepsilon_{i} \\ \end{vmatrix} = 0 \] 可得 \[ \left \{ \begin{aligned} & \varepsilon_{1} = -0.918 \\ & \varepsilon_{2} = 2.810 \\ \end{aligned} \right. \] 将能量较低的\(\varepsilon_{1}\)代入 \[ c_{11}(F_{21} - \varepsilon_{1}S_{21}) + c_{21}(F_{22} - \varepsilon_{1}S_{22}) = 0 \] 可得 \[ \dfrac{c_{11}}{c_{21}} = 4.62 \] 由归一化不难得到 \[ c_{11}^{2} + c_{21}^{2} + 2c_{11}c_{21}S_{12} = 1 \]\(S_{12} = 0.8366\)以及\(\dfrac{c_{11}}{c_{21}} = 4.62\)代入可得 \[ \left \{ \begin{aligned} & c_{11} = 0.842 \\ & c_{21} = 0.182 \\ \end{aligned} \right. \]

第五次计算密度矩阵

新密度矩阵的矩阵元分别为 \[ \left \{ \begin{aligned} & P_{11} = 2c_{11}^{2} = 1.418 \\ & P_{12} = P_{21} = 2c_{11}c_{21} = 0.306 \\ & P_{22} = 2c_{21}^{2} = 0.066 \\ \end{aligned} \right. \]

不难看出,在当前精度下,我们可以认为此时SCF已经收敛。由此可得计算能量

能量计算

由于Hartree-Fock能量可以表示为 \[ E_{\mathrm{HF}} = \dfrac{1}{2}\sum\limits_{r = 1}^{b}\sum\limits_{s = 1}^{b}P_{rs}(F_{rs} + H_{rs}^{\mathrm{core}}) + V_{NN} \] 对于He原子而言,\(V_{NN} = 0\)。将\(b = 2\)和其它数据代入可得 \[ E_{\mathrm{HF}} = -2.862\ \mathrm{Hartree} \] 此时\(\mathrm{Hartree}\)为原子单位下的能量单位,其大小为\(4.3598 \times 10^{-18}\ \mathrm{J}\)。由此 \[ E_{\mathrm{HF}} = -1.248\times 10^{-17}\ \mathrm{J} = -77.89\ \mathrm{eV} \]

参考文献

[1] Ira N. Levine.Quantum Chemistry[M].Seventh edition.Boston:Pearson,2014:407-416.

[2] Attila Szabo, Neil S. Ostlund.Modern Quantum Chemistry: introduction to advanced electronic structure theory[M].First edition.Mineola, N.Y.:Dover Publications,1996:142-145.

[3] Henry Eyring, The Late John Walter, George E. Kimball.Quantum Chemistry[M].First edition.John Wiley & Sons, Inc.,1944:369-371.